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1  SUMMARY 


Terahertz  (THz)-band  eommunieation  is  an  emerging  teehnology  that  ean  address  two  relevant 
ehallenges  for  the  U.S.  Air  Foree.  First,  THz-band  frequeneies  (0.1  to  lOTHz)  remain  almost 
eompletely  unutilized  for  eommunieation;  henee,  being  the  first  to  exploit  this  band  provides  an 
unmatehed  teehnologieal  advantage.  Seeond,  the  huge  transmission  bandwidth  assoeiated  with  the 
THz-band  ean  satisfy  demand  for  higher  data  rates  in  eommunieation  networks,  thereby  enabling 
Terabit-per- seeond  (Tbps)  wireless  links.  The  objeetive  of  this  projeet  is  to  establish  the  theoretieal 
foundations  for  praetieal  THz-band  eommunieation  networks.  Over  the  duration  of  the  projeet  (20 
months),  eontributions  have  been  made  aeross  different  areas  both  related  to  THz  deviee  teehnologies 
as  well  as  to  physieal  and  link  layer  solutions  for  THz  eommunieation  networks.  More  speeifieally, 
in  terms  of  THz  deviees,  graphene-based  plasmonie  phase  modulators  and  plasmonie  nano-antenna 
arrays  have  been  proposed,  modeled  and  analyzed.  In  terms  of  physieal  and  link  layer  solutions,  new 
models  for  interferenee  in  THz  eommunieation,  symbol  deteetion  and  physieal  layer  synehronization 
strategies;  a  paeket  size  optimization  framework;  and  link  layer  synehronization  and  medium  aeeess 
eontrol  protoeols  have  been  devised,  formulated  and  modeled.  These  have  also  been  validated  by  a 
eombination  of  numerieal,  simulation  and  experimental  results. 
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2  INTRODUCTION 


Over  the  last  few  years,  wireless  data  traffie  has  drastieally  inereased  due  to  a  ehange  in  the  way 
today’s  soeiety  ereates,  shares  and  eonsumes  information.  In  2015,  the  need  to  provide  wireless 
eonneetivity  to  anything,  anywhere,  anytime  has  resulted  in  more  than  7.9  million  mobile  deviees 
eonneeted  to  the  Internet,  whieh  have  generated  a  total  of  3.7  exabytes  per  month  of  mobile  data 
traffie  [1].  Moreover,  estimates  foreeast  that  there  will  be  11.6  billion  mobile-eonneeted  deviees  by 
2020,  mainly  in  part  due  to  the  Internet  of  Things  paradigm.  In  parallel  with  the  growth  in  the  total 
number  of  intereonneeted  deviees,  there  has  been  an  inereasing  demand  for  higher  speed  wireless 
eommunieation.  In  partieular,  wireless  data  rates  have  doubled  every  eighteen  months  for  the  last 
three  deeades  [2].  Following  this  trend,  Terabit-per-second  (Tbps)  links  are  expeeted  to  beeome  a 
reality  within  the  next  hve  years. 

Several  alternatives  are  being  eonsidered  to  meet  this  demand.  At  frequencies  below  5  GHz, 
advaneed  digital  modulations  and  sophistieated  physieal  layer  sehemes  are  being  used  to  aehieve 
a  very  high  speetral  effieieney.  However,  the  very  small  available  bandwidth  in  the  overerowded 
eleetromagnetie  (EM)  speetrum  limits  the  aehievable  data  rates.  For  example,  in  Long-Term 
Evolution  Advaneed  (LTE-A)  networks,  peak  data  rates  on  the  order  of  1  Gigabit-per-seeond 
(Gbps)  are  feasible  over  100  MHz  of  aggregated  bandwidth  when  utilizing  a  4x4  MIMO  seheme  [3]. 
Similarly,  dynamie  speetrum  aeeess  and  sharing  teehniques  are  being  heavily  investigated  to  better 
utilize  the  available  bandwidth.  While  these  teehnologies  will  eontribute  to  aehieve  maximum 
speetral  effieieney,  they  are  far  from  being  able  to  support  Tbps. 

Millimeter-wave  (mm-wave)  communication  systems  (30  to  300  GHz)  have  gained  a  lot  of 
attention  in  the  last  few  years  due  to  their  ability  to  support  mueh  higher  data  rates  than  eommunieation 
systems  below  5  GHz  [4].  Several  sub-bands  have  been  alloeated  for  eommunieations,  ineluding 
38.6  to  40  GHz  (1.4  GHz  bandwidth),  57  to  64  GHz  (7  GHz  bandwidth,  but  usually  smaller  due  to 
international  regulations)  and  71  to  76  GHz  /  8 1  to  86  GHz  (10  GHz  bandwidth  in  total).  Millimeter- 
wave  teehnology  is  already  playing  a  key  role  in  eurrent  systems  (e.g.,  IEEE  802.1  lad),  and  will 
eontinue  to  do  so  in  future  (e.g.,  5G  eellular)  systems.  While  this  is  the  way  to  go,  the  total 
eonseeutive  available  bandwidth  for  mm-wave  eommunieation  systems  is  still  less  than  10  GHz. 
This  would  require  a  physieal  layer  effieieney  of  almost  100  bps/Hz  to  support  Tbps,  whieh  is  several 
times  higher  than  that  of  the  state  of  the  art  for  existing  systems.  This  result  motivates  further 
researeh  and  the  exploration  of  even  higher  frequeney  bands  for  eommunieation. 

Optical  wireless  communication  ( OWC)  systems,  whieh  operate  in  the  infrared  (187  to  400  THz 
/  750  to  1600  nm),  visible  (400  to  770  THz  /  390  to  750  nm),  or  even  ultraviolet  (1000  to  1500  THz  / 
200  to  280  nm)  EM  speetrum  bands,  are  similarly  being  explored  as  a  way  to  improve  the  aehievable 
data  rates  in  wireless  networks  [5].  The  intrinsieally  very  large  available  bandwidth  at  sueh  very 
high  frequeney  plays  to  their  advantage.  However,  there  are  several  aspeets  that  eurrently  limit 
these  approaehes  and  require  further  researeh,  ineluding  the  size  and  limited  portability  of  infrared 
systems;  low  transmission  power  budget  due  to  eye-safety  limits  and  the  impaet  of  atmospherie 
effeets  on  the  signal  propagation  (e.g.,  seattering  from  moleeules  and  partieles);  or  the  impaet  of 
ambient  noise  in  visible  eommunieation  (e.g.,  any  light  souree). 

In  this  eontext.  Terahertz-band  (0.1  to  10  THz)  communication  is  envisioned  as  a  key 
wireless  technology  to  satisfy  the  need  for  much  higher  wireless  data  rates  [6-9].  This 
frequency  band,  which  lies  between  mm-waves  and  the  far  infrared,  remains  still  one  of  the  least 
explored  regions  in  the  EM  spectrum.  For  many  decades,  the  lack  of  compact  high-power  signal 
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sources  and  high- sensitivity  detectors  able  to  work  at  room  temperature  has  hampered  the  use  of  the 
THz  band  for  any  application  beyond  sensing.  However,  many  recent  advancements  with  different 
technologies  is  finally  closing  the  so-called  THz  gap.  For  example,  on  the  one  hand,  in  an  electronic 
approach,  TTT-V  semiconductor  technologies  have  demonstrated  record  performance  in  terms  of 
output  power,  noise  figure,  and  power-added  efficiency  at  sub-THz  frequencies,  and  are  quickly 
approaching  the  1  THz  mark  [10-12].  On  the  other  hand,  in  an  optics  approach.  Quantum  Cascade 
Lasers  are  rising  as  potential  candidates  for  high-power  THz -band  signal  generation  [13-16].  More 
recently,  the  use  of  nanomaterials  such  as  graphene  is  enabling  the  development  of  novel  plasmonic 
devices,  which  intrinsically  operate  in  the  THz-band  [17-22]. 

Independently  of  the  specific  enabling  technology,  THz-band  communication  can  address 
two  relevant  challenges  for  the  US  Air  Force  (USAF).  First,  THz-band  frequencies  remain  almost 
completely  unutilized  for  communication;  hence,  being  the  first  to  exploit  this  band  provides  an 
unmatched  technological  advantage.  Second,  the  huge  transmission  bandwidth  associated  with  the 
THz-band  can  satisfy  demand  for  higher  data  rates  in  communication  networks,  thereby  enabling 
Ter abit-per- second  wireless  links  [8].  Nevertheless,  this  very  large  bandwidth  comes  at  the  cost  of  a 
very  high  propagation  loss  [23],  mainly  because  of  molecular  absorption,  which  also  creates  a  unique 
distance  dependence  on  the  available  bandwidth.  All  these  introduce  many  challenges  to  practical 
THz  communication  systems  and  require  the  development  of  innovative  solutions.  Moreover,  many 
of  these  might  be  helpful  for  broadband  wireless  communication  systems  below  and  above  the  THz 
band,  i.e.,  mm-waves  and  OWC,  respectively. 

The  objective  of  this  project  is  to  establish  the  theoretical  foundations  for  practical  THz-band 
communication  networks.  Along  this  20-month-long  project,  contributions  have  been  made  across 
different  areas,  ranging  from  new  graphene-based  plasmonic  devices  for  true  THz  communications, 
to  new  physical  and  link  layer  protocols  for  ultra-broadband  networking  in  the  THz  band.  More 
specifically,  contributions  have  been  made  along  the  following  intertwined  thrusts: 

•  Graphene-based  Plasmonic  Phase  Modulator  (Sec.  3):  While  major  efforts  have  been 
made  towards  developing  compact  THz  sources  and  detectors  able  to  work  at  room  temperature, 
the  modulator/demodulator  is  a  fundamental  block  in  any  communication  system.  As  part  of 
this  effort,  we  have  proposed,  modeled  and  analyzed  the  performance  of  a  graphene -based 
modulator  able  to  tune  on  demand  the  phase  of  a  propagating  Surface  Plasmon  Polariton  (SPP) 
wave  as  it  propagates  from  the  source  towards  the  antenna.  This  work  has  been  published 
in  [24]  and  an  extended  version  of  the  work  is  being  prepared  for  journal  publication.  In 
addition,  the  University  at  Buffalo  is  currently  preparing  a  US  Patent  Application  for  this 
invention. 

•  Graphene-based  Plasmonic  Nano-antenna  Arrays  (Sec.  4):  The  limited  power  of  THz 
transceivers  and  the  very  high  path  loss  at  THz  frequencies  drastically  limit  the  transmission 
distance  between  THz  devices.  To  overcome  this  limitation,  by  taking  advantage  of  the 
very  small  size  of  THz  plasmonic  transceivers  and  antennas,  massive  nano-antenna  arrays 
can  be  created  [25].  As  part  of  this  project,  we  have  proposed,  modeled  and  analyzed  the 
performance  of  graphene -based  plasmonic  nano-antenna  arrays,  whose  radiation  diagram 
can  be  dynamically  tuned  in  space  and  frequency.  This  work  has  been  published  in  [26], 
and  an  extended  version  of  the  work  is  being  prepared  for  journal  publication.  In  addition, 
new  strategies  to  experimentally  fabricate  the  proposed  structures  have  been  devised  and 
disseminated  to  the  research  community  [27].  Moreover,  we  have  set  up  a  pioneering 
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experimental  eharaeterization  platform  for  2D  nanomaterials  and  nanostruetures,  to  be  used 
for  the  analysis  of  the  developed  deviees,  and  published  the  results  in  [28]. 

Interference  Modeling  in  THz  Pulse-based  Communications  (Sec.  5):  Many  of  the 
envisioned  applieations  of  THz -band  eommunieations  involved  a  large  number  of  deviees 
eommunieating  over  short  distanees  at  very  high  data-rates  (up  to  Tbps)  [29].  In  sueh  seenarios, 
multi-user  interferenee  ean  beeome  a  bottleneek  if  not  properly  addressed.  As  part  of  this 
effort,  we  have  developed  and  experimentally  validated  a  stoehastie  multi-user  interferenee 
model  for  THz-band  eommunieation  networks.  The  proposed  model  takes  into  aeeount  not 
just  the  power  of  interfering  signals,  but  also  their  preeise  time-domain  waveforms.  An 
innovative  testbed  has  been  developed  to  validate  the  analytieal  models.  This  work  has  been 
published  in  [30]  and  an  extended  version  is  eurrently  in  preparation.  In  addition,  we  have 
also  explored  and  published  both  stoehastie  geometry  [31]  and  novel  integral  geometry  [32] 
approaehes  to  interferenee  modeling  in  THz  eommunieation  networks. 

Symbol  Detection  and  Physical  Layer  Synchronization  (Sec.  6):  When  transmitting  at 
Gbps  and  Tbps,  physieal  layer  synehronization  beeomes  one  of  the  main  ehallenges  due  to  the 
very  short  symbol  duration.  Due  to  the  laek  of  Analog-to-Digital  Converters  (ADCs)  able  to 
digitize  the  reeeived  signal  at  rates  elose  to  the  Nyqvist  frequeney,  analog  synehronization 
sehemes  are  needed.  As  part  of  this  effort,  we  have  proposed,  modeled  and  experimentally 
tested  an  analog  symbol  deteetion  and  synehronization  seheme,  built  upon  graphene -based 
voltage-eontrolled  delay  (VCD)  lines  and  a  eontinuous-time  moving- average  (CTMA)  deteetor. 
The  developed  worked  has  been  published  in  [33]  and  an  extended  version  is  in  preparation. 

Packet  Size  Optimization  for  THz  Communication  (Sec.  7):  One  of  the  key  applieations 
of  THz  eommunieations  is  in  the  field  of  nanonetworks  and  wireless  nanosensor  networks 
(WNSNs)  [29,34].  While  the  THz  band  provides  nanosensors  with  a  very  large  bandwidth  [23], 
the  limited  eapaeity  of  nano-batteries,  whieh  requires  the  use  of  time-eonsuming  energy¬ 
harvesting  proeedures  [35,36],  and  the  limited  eomputational  eapabilities  of  nanosensors, 
affeet  the  throughput  of  WNSNs.  All  these  interdependeneies  motivate  the  joint  analysis  of 
the  nano-deviee  eapabilities,  the  THz  band  peeuliarities  and  their  impaet  on  the  aehievable 
throughput.  As  part  of  this  effort,  we  have  taekled  the  throughput  maximization  problem  in 
WNSNs,  by  taking  into  aeeount  the  existing  deviee  and  eommunieation  interdependeneies, 
and  identified  the  optimal  paeket  size  for  different  system-level  objeetives.  The  developed 
work  has  been  published  in  [37]. 

Link  Layer  Synchronization  and  Medium  Access  Control  Protocols  for  THz  Networks 
(Sec.  8):  Traditionally,  Medium  Aeeess  Control  (MAC)  protoeols  have  been  designed  by 
taking  into  aeeount  that  the  most  searee  resouree  is  the  ehannel  available  bandwidth.  As  a 
result,  nodes  usually  need  to  eontend  for  the  ehannel.  At  THz  frequeneies,  perhaps  the  only 
resouree  whieh  is  not  limiting  is  the  bandwidth.  However,  as  mentioned  earlier,  this  bandwidth 
eomes  at  the  eost  of  a  very  high  propagation  loss,  whieh  imposes  both  energy  eonstraints  in 
battery-powered  deviees  and  the  need  for  ultra-direetional  eommunieations  among  spaeed 
deviees.  As  part  of  this  effort,  we  have  proposed,  modeled  and  analyzed  the  performanee  of 
new  MAC  protoeols  that  shift  the  foeus  from  the  ehannel  to  the  reeeiver.  The  proposed  reeeiver- 
initiated  protoeols  rely  on  a  one-way  handshake  triggered  by  the  reeeiver,  and  ean  effeetively 


Approved  for  Public  Release;  Distribution  Unlimited 

4 


outperform  existing  zero- way  (i.e.,  ALOHA-type)  or  2-way  (i.e.,  CSMA-type)  protoeols  in 
terms  of  aehievable  throughput.  This  work  has  been  published  in  [38],  and  an  extended 
journal  version  has  been  submitted  for  publieation.  In  addition,  we  have  proposed,  modeled 
and  analyzed  a  dual-band  MAC  protoeol,  aimed  at  synergistieally  leveraging  the  advantages 
of  both  2.4  GHz  (omnidireetional  transmissions)  and  THz  (ultra-high  data-rates)  [39]. 

In  the  following  seetions,  we  deseribe  the  Methods,  Assumptions,  Proeedures,  Results  and 
Diseussions  for  eaeh  of  the  aforementioned  eontributions. 
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3  GRAPHENE-BASED  PLASMONIC  PHASE  MODULATOR 


3.1  Motivation 

In  the  last  decade,  there  has  been  major  progress  towards  developing  new  types  of  THz  sources 
and  detectors  [10-16,21],  usually  in  the  context  of  THz  sensing  and  imaging  applications.  When  it 
comes  to  communications,  in  addition  to  THz  signal  sources  and  detectors,  a  modulator  is  needed  to 
embed  information  on  the  transmitted  signals.  The  desired  properties  of  a  modulator  include  high 
modulation  bandwidth,  i.e.,  the  speed  at  which  the  properties  of  the  modulated  signal  can  be  changed, 
and  high  modulation  depth,  i.e.,  the  maximum  difference  between  modulation  states.  Different 
types  of  modulators  able  to  control  the  amplitude  or  phase  of  THz  waves  have  been  developed  to 
date  [40].  In  [41],  a  high-electron-mobility  transistor  based  on  a  III-V  semiconductor  material  was 
utilized  to  modulate  the  amplitude  of  a  THz  wave.  In  [42],  a  metamaterial-based  modulator  was 
utilized  to  control  the  phase  of  a  THz  wave.  In  both  cases,  sub-GHz  modulation  bandwidths  and 
low  modulation  depths  limit  the  use  of  these  devices  in  practical  communication  systems. 

More  recently,  the  use  of  graphene  to  develop  THz  wave  modulators  has  been  proposed  [43]. 
Graphene  has  excellent  electrical  conductivity,  making  it  very  well  suited  for  propagating  extremely- 
high-frequency  electrical  signals  [18].  In  [44],  a  graphene-based  amplitude  modulator  for  THz 
waves  was  developed.  This  was  enabled  by  the  possibility  to  dynamically  control  the  conductivity 
of  graphene.  In  [45],  a  similar  principle  was  exploited  in  a  graphene-based  meta-device.  In  these 
setups,  the  main  challenge  is  to  increase  the  modulation  depth.  A  low  modulation  depth  makes 
the  transmitted  symbols  more  difficult  to  distinguish  and,  thus,  results  in  higher  symbol  error  rates 
(SER)  in  practical  communication  systems. 

As  part  of  this  project,  we  have  proposed,  modeled  and  analyzed  the  performance  of  a  graphene- 
based  plasmonic  phase  modulator  for  THz -band  communications.  Our  proposed  modulator  consists 
of  a  fixed-length  graphene -based  plasmonic  waveguide  with  a  metallic  back  gate.  Its  working 
principle  relies  on  the  possibility  to  electronically  control  the  propagation  speed  of  a  Surface  Plasmon 
Polariton  (SPP)  wave  on  graphene  at  THz  frequencies  by  modifying  the  chemical  potential  of  the 
graphene  layer.  In  this  section,  we  summarize  our  contributions  towards  this  aim.  In  Sec.  3.2.1, 
starting  from  the  dynamic  complex  conductivity  of  graphene  and  a  revised  dispersion  equation 
for  SPP  waves,  we  develop  an  analytical  model  for  the  plasmonic  phase  modulator.  By  utilizing 
the  model,  in  Sec.  3.2.3,  we  analyze  the  performance  of  the  proposed  plasmonic  modulator  when 
utilized  to  implement  a  M-ary  phase  shift  keying  modulation  in  terms  of  SER.  Einally,  we  validate 
the  model  by  means  of  electromagnetic  simulations,  and  provide  numerical  results  to  illustrate  the 
modulator  performance  in  Sec.  3.3. 

3.2  Methods,  Assumptions  and  Procedures 

3.2.1  Propagation  Properties  of  SPP  Waves  on  Gated  Graphene  Structures 

The  analysis  of  the  proposed  plasmonic  phase  modulator  requires  the  characterization  of  the  SPP 
propagation  properties  on  graphene.  These  depend  on  the  conductivity  of  the  graphene  sheet.  In  this 
section,  we  first  recall  the  conductivity  model  utilized  in  our  analysis  and  then  define  the  dispersion 
equation  for  SPP  waves  on  gated  graphene  structures. 
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Complex  Conductivity  Model  of  Graphene  In  our  analysis,  we  consider  a  surface  conductivity 
model  for  infinitely  large  graphene  sheets  obtained  using  the  Kubo  formalism  [46,47].  This  is  given 
by 
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where  oj  =  2nf,  fi  =  hjlTi  is  the  reduced  Planck’s  constant,  e  is  the  electron  charge,  kB  is  the 
Boltzmann  constant,  T  is  temperature,  Ep  refers  to  the  Fermi  energy  of  the  graphene  sheet,  and 
Tg  is  the  relaxation  time  of  electrons  in  graphene,  which  depends  on  the  electron  mobility  Hg.  Ep 
can  be  easily  modified  by  means  of  electrostatic  bias  or  gating  of  the  graphene  layer,  enabling  the 
aforementioned  antenna  tuning. 

As  we  showed  in  [19],  a  more  accurate  conductivity  model  can  be  developed  by  taking  into 
account  the  impact  of  electron  lateral  confinement  on  graphene  nano-ribbons,  but  the  two  models 
converge  for  graphene  strips  which  are  50  nm  wide  or  more.  In  our  analysis,  we  consider  plasmonic 
resonant  cavities  which  are  a  few  hundred  nanometers  wide.  Finally,  from  [46,47],  it  is  important  to 
note  that  the  conductivity  model  described  by  (1)  and  the  following  was  derived  by  neglecting  the 
spatial  dispersion  of  the  AC  field.  Therefore,  it  can  be  used  for  the  analysis  of  the  SPP  propagation 
in  the  long  wavelength  limit  only,  i.e.,  co  »  ksppVp,  where  kspp  is  the  SPP  wave  number  and 
y/r  «  8  X  10^  m/s  is  the  Fermi  velocity  of  Dirac  fermions  in  graphene. 

Dispersion  Equation  for  Surface  Plasmon  Polariton  Waves  The  propagation  properties  of  SPP 
waves  can  be  obtained  by  deriving  and  solving  the  SPP  wave  dispersion  equation  on  graphene.  In 
many  of  the  related  graphene  plasmonic  works  [48-50],  the  dispersion  equation  was  obtained  by 
considering  a  graphene  layer  at  the  interface  between  two  infinitely  large  dielectric  materials,  usually 
between  air  and  silicon  dioxide  (Si02).  However,  the  proposed  modulator  relies  on  the  presence 
of  a  metallic  ground  plane  at  a  distance  d  from  the  graphene  layer,  which  is  needed  both  to  create 
the  plasmonic  waveguide  as  well  as  to  control  the  Fermi  energy  of  the  graphene  layer  and  tune  its 
conductivity. 

From  [51],  the  dispersion  equation  for  Transverse  Magnetic  (TM)  SPP  waves  on  gated  graphene 
structures  in  the  quasi-static  regime — i.e.,  for  kspp  »  m/c,  where  c  is  the  speed  of  light — is  given  by 


_  cr^  _  SI+S2  coth  {ksppd) 
ojsq  kspp 


(5) 


where  cr^  is  the  conductivity  of  graphene  given  by  (1),  s\  is  the  relative  permittivity  of  the  dielectric 
above  the  graphene  layer,  and  S2  is  the  relative  permittivity  of  the  dielectric  between  the  graphene 
layer  and  the  metallic  ground  plane,  which  are  separated  by  a  distance  d.  It  can  be  easily  shown  by 
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taking  the  limit  of  J  ^  oo  that  (5)  tends  to  the  quasi-statie  dispersion  equation  of  SPP  waves  in 
ungated  graphene  used  in  the  aforementioned  works. 

By  solving  (5),  the  eomplex  wave  veetor  kspp  ean  be  obtained.  The  real  part  of  the  wave 
veetor, 
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determines  the  SPP  wavelength  Aspp  and  the  SPP  wave  propagation  speed.  The  imaginary  part 
determines  the  SPP  dee  ay  or,  inversely, 


llmiksppY 

determines  the  SPP  propagation  length,  whieh  is  dehned  as  the  distanee  at  whieh  the  SPP  intensity 
has  deereased  by  a  value  of  1/e.  Unfortunately,  a  elosed-form  expression  for  kspp  in  this  ease  does 
not  exist,  but  ean  only  be  obtained  numerieally. 

3.2.2  Graphene-based  Plasmonic  Phase  Modulator 

In  this  seetion,  we  explain  the  working  prineiple  of  the  plasmonie  phase  modulator  and  develop  its 
analytieal  model. 

Working  Principle  The  eoneeptual  design  of  the  proposed  graphene-based  plasmonie  phase 
modulator  is  shown  in  Fig.  1.  The  phase  modulator  eonsists  of  a  plasmonie  waveguide,  whieh  is 
eomposed  of  a  graphene  sheet  (the  plasmonie  material)  mounted  over  a  metallie  flat  surfaee  (the 
ground  plane),  with  a  dieleetrie  material  layer  in  between,  whieh  supports  the  graphene  layer.  In 
this  paper,  we  eonsider  that  a  SPP  wave  is  already  propagating  over  the  graphene  layer.  Different 
meehanisms  eould  be  utilized  to  launeh  the  SPP  wave.  For  example,  a  THz  plasma  wave  eould  be 
generated  by  a  HEMT-like  deviee  and  eoupled  to  the  graphene  waveguide  [21].  Alternatively,  a 
Quantum  Caseade  Laser  and  a  grating  strueture  eould  be  utilized  for  the  same  [14]. 


Fig.  1 :  Graphene-based  plasmonic  phase  modulator. 


The  basie  idea  for  any  phase  modulator  is  to  establish  a  relation  between  the  data  bits  to 
be  transmitted,  whieh  eonstitute  the  modulating  signal,  and  the  transmitted  signal  phase.  In  our 
proposed  plasmonie  phase  modulator,  the  modulating  signal  is  applied  as  a  bias  voltage  to  the 
graphene  layer  and  eontrols  its  Fermi  energy,  Ep.  From  (6),  (5)  and  (1),  it  is  elear  that  the  Fermi 
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energy  direetly  eontrols  the  propagation  speed  of  the  SPP  wave  on  the  gated  graphene  strueture. 
Therefore,  for  a  fixed  length  waveguide,  the  phase  at  the  output  of  the  waveguide  is  effeetively 
modulated  by  the  data  bits  to  be  transmitted.  The  resulting  signal  ean  be  further  propagated  and 
eventually  radiated  in  free-spaee  by  means  of  a  plasmonie  nano-antenna  [19,49, 50]. 

Analytical  Model  In  this  seetion  we  develop  an  analytieal  model  for  the  phase  of  the  plasmonie 
signal  at  the  output  of  the  plasmonie  phase  modulator  shown  in  Fig.  1.  We  denote  the  plasmonie 
signals  at  the  input  and  the  output  of  the  plasmonie  waveguide  as  X  and  T,  respeetively.  The 
modulator  frequeney  response  is  denoted  by  H.  The  following  relation  ean  be  then  written, 

Y(f,EF)  =  X(f)H(f,EF),  (8) 

where  /  stands  for  frequeney  and  Ef  is  the  Fermi  energy  of  the  graphene  layer  on  whieh  the  SPP 
wave  propagates. 

The  modulator  frequeney  response  H  is  given  by 

H  (/,  Ef)  =  \H  (/,  Ef)\  exp  (j9  (/,  Ef)),  (9) 

where  \H\  aeeounts  for  the  variation  in  the  SPP  wave  intensity  and  6  represents  the  ehange  in  the 
SPP  phase  at  the  output  of  the  fixed-length  waveguide. 

From  See.  3.2.1,  the  magnitude  of  the  modulator  response  ean  be  written  as 

\H(f,EF)\  =  exp  {-2lm{kspp(f,EF)}L),  (10) 

where  L  represents  the  waveguide  length. 

The  total  phase  ehange  9  that  the  SPP  wave  suffers  as  it  propagates  through  the  waveguide  is 
given  by 

'2,71  L 

^ =  ,  rr-Fn  =  ^ (/> Ef)  ,  (11) 

'ispp  [j  ,  Ej) 

where  Aspp  is  the  plasmonie  wavelength  obtained  from  kspp  as  diseussed  in  See.  3.2.1,  whieh  depends 
on  the  signal  frequeney  /  and  the  Fermi  energy  Ef- 

By  eombining  (10)  and  (1 1)  in  (9),  the  modulator  frequeney  response  ean  be  written  as 

H(f,EF)  =cxp{-2lm{k,pp(f,EF)}L) 

•  exp  (j  Re  {k,pp  (/,  Ef)}  L)  . 

In  an  ideal  phase  modulator,  the  intensity  or  amplitude  of  the  signal  should  remain  eonstant, 
independently  of  the  phase.  However,  the  SPP  deeay  in  graphene  struetures  is  not  negligible.  As  a 
result,  we  eannot  independently  modulate  the  signal  amplitude  and  phase.  This  has  a  direet  impaet 
on  the  performanee  of  the  modulator  in  a  praetieal  eommunieation  system,  whieh  we  analyze  in  the 
next  seetion. 

3.2.3  Performance  Analysis 

In  this  section,  we  define  the  constellation  of  a  non-uniform  plasmonie  phase  shift  keying  digital 
modulation  and  formulate  the  SER  for  M-ary  modulations. 
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Signal  Space  Constellation  The  signal  space  or  constellation  represents  the  possible  symbols 
generated  by  a  given  modulation  scheme  as  points  in  the  complex  plane.  The  real  part  of  each  of 
such  points  is  referred  to  as  the  in-phase  component  and  the  imaginary  part  denotes  the  quadrature 
component. 

The  number  of  modulated  symbols  or  points  in  the  constellation  is  given  by  M  =  2^,  where 
k  =  2,  A, .. .  refers  to  the  modulation  order.  The  position  of  each  symbol  Sm,  m  =  1..M,  depends  on 
the  modulator  behavior.  For  the  system  described  in  Sec.  3.2.2,  at  fixed  carrier  frequency  /c,  the 
magnitude  and  phase  of  each  symbol  is  given  by 

>5™  =  |S„|exp(0„),  (13) 

\Sm\=AQ\H{f„Ep^m)l  (14) 

0m  =  Oo  +  6  [fc,  Ef^m)  ,  (15) 

where  Ao  and  00  refer  to  the  amplitude  and  phase  of  the  input  SPP  wave.  Ep^m  = 
is  the  set  of  Fermi  energies  that  correspond  to  the  transmitted  symbols.  In  our  analysis,  we  consider 
Ao  =  1  and  6q  =  0.  The  constellation  for  the  proposed  plasmonic  phase  modulator  can  only  be 
numerically  obtained  and  will  be  provided  in  Sec.  3.3. 

Symbol  Error  Rate  The  most  common  metric  for  a  modulation  scheme  in  a  practical  communi¬ 
cation  system  is  the  SER.  This  is  implicitly  related  to  the  modulation  intensity  or  depth.  The  more 
“distinguishable"  the  symbols  are,  the  lower  the  SER.  In  general  terms,  for  a  modulated  symbol  Sm, 
the  symbol  error  probability  Pg  is  given  by  [52], 

Pg  =  P  {Detect  Sm,m  m  \  Given  that  Sm  is  sent} ,  (16) 


where  m  =  1, 2, 3, . . . ,  M.  The  SER  for  a  digital  phase  modulation  with  uniform  constellation  is 
derived  based  on  the  symbol  decision  regions,  which  due  to  symmetry,  are  easy  to  define.  However, 
this  is  not  the  case  for  non-uniform  modulations.  Instead,  the  SER  for  the  proposed  plasmonic 
phase  modulation  scheme  can  be  directly  derived  starting  from  the  distance  between  symbols  in  the 
non-uniform  constellation. 

In  general  terms,  the  union  bound  for  the  SER  is  given  by  [52], 


SER  < 


1 

M 


M 

E  E  fi 

m=l  \<m<M 
ih+m 


I  D{Sm,  Sfhf' 

~Wo 


(17) 


where  the  Q  function  refers  to  the  tail  probability  of  the  standard  normal  distribution,  D{Sm,  Sm) 
stands  for  the  distance  between  two  symbols  Sm  and  Sm,  and  is  given  by 


D{Sm,  Sm)  —  \\Sm  Sm\ 


(18) 


and  Nq  is  the  noise  power  spectral  density. 

A  common  representation  of  the  SER  is  as  a  function  of  signal-to-noise  ratio  (SNR)  or  the 
energy  per  symbol  to  noise  power  spectral  density  Es/Nq.  Erom  (14),  this  is  given  by 


SNRm 


No 


Hm{fc,EF)\^ 

No 


(19) 
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Finally,  by  combining  (17),  (18)  and  (19),  the  SER  for  the  non-uniform  eonstellation  ean  be 
further  expressed  as 


m=l  l<m<M 
ih+m 

The  SER  will  be  numerieally  investigated  in  the  next  seetion. 

3.3  Results  and  Discussions 

In  this  seetion,  we  validate  our  models  and  analyze  the  performanee  of  the  proposed  plasmonie 
phase  modulator. 

3.3.1  Model  Validation 

We  utilize  COMSOL  Multi-physies  to  simulate  the  behavior  of  the  plasmonie  phase  modulator 
shown  in  Eig.  1 .  Graphene  is  modeled  as  a  transition  boundary  eondition  with  impedanee  given 
by  (1),  with  Tg  =  2.2  ps  at  room  temperature  T  =  300  K.  The  graphene  layer  rests  on  top  of  a  metallie 
ground  plane  with  a  90  nm-thiek  Si02  dieleetrie  in  between  {sr  =  4).  In  Eig.  2,  the  z-eomponent 
of  the  eleetrie  field  on  a  2-/im-long  graphene-based  waveguide  is  shown  for  /c=4  THz  and  two 
values  of  Ep,  namely,  0.13  eV  and  0.28  eV.  Eor  Ep  =  0.13  eV,  the  waveguide  length  L  eorresponds 
to  approximately  one  and  a  half  SPP  wavelength  and,  thus,  it  introduees  a  phase  ehange  of  n.  Eor 
Ep  =  0.28  eV,  the  waveguide  length  L  eorresponds  to  one  full  SPP  wavelength  and  introduees  a 
phase  ehange  of  2n.  Henee,  we  ean  define  a  plasmonie  phase  modulator  of  order  M=2,  where  bit 
“0"  is  transmitted  as  a  phase  ehange  of  In  by  tuning  Ep  to  0.28  eV  and  bit  “1"  is  transmitted  with  a 
phase  ehange  of  n  by  tuning  to  0.13  eV. 

In  addition  to  the  phase,  we  are  interested  in  the  ehange  in  the  amplitude  of  the  SPP  wave 
amplitude,  as  it  will  affeet  the  signal  spaee  eonstellation  and  the  SER.  In  Eig.  3(a)  and  Eig.  3(b),  we 
illustrate  the  magnitude  and  phase  of  the  plasmonie  phase  modulator  as  a  funetion  of  the  Eermi 
energy,  Ep,  and  for  different  earner  frequeneies,  fc.  On  the  one  hand,  we  are  interested  in  working 
in  a  range  of  Ep  sueh  that  the  magnitude  of  the  modulator  does  not  signifieantly  ehange.  On  the 
other  hand,  however,  we  need  at  least  a  phase  differenee  of  n  to  ereate  orthogonal  symbols.  Next, 
we  investigate  the  performanee  of  a  speeifie  modulator  design. 

3.3.2  Constellation  and  Symbol  Error  Rate 

In  Eig.  4(a),  the  non-uniform  eonstellation  for  a  plasmonie  phase  modulator  with  /c  =  4  THz, 
Epfi  =  0.28  eV,  Ep^i  =  0.13eV,  and  L  =  2  pm,  d  =  90  nm  is  shown.  Similarly,  in  Eig.  4(b),  the 
SER  (17)  for  the  proposed  modulator  is  shown  as  a  funetion  of  the  SNR  and  eompared  to  that  of  a 
uniform  binary  phase  shift  keying  (BPSK)  modulation  with  the  same  average  energy  per  symbol  Es. 
As  expeeted,  the  SER  for  the  proposed  modulator  is  slightly  higher  but  still  eomparable  to  that  of 
the  uniform  ease. 
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(a)  Ep  -  0.13  eV 
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(b)  Ep  =  0.28  eV 


Fig.  2:  Electric  field  distribution  over  a  graphene-based  waveguide  at  /c=4THz,  for  different  Fermi  energies, 
Ep  (L=2pm,  d=90nm). 


(a)  Magnitude,  \H{fc,Ep)\ 


(b)  Phase,  6{fc,Ep) 


Fig.  3:  Frequency  response  of  the  phase  modulator,  as  a  function  of  the  Fermi  energy,  Ep,  for  different  carrier 
frequencies  fc  (L=2  pm,  r/=90nm). 

3.3.3  Discussion 


As  part  of  this  effort,  we  have  proposed  a  novel  graphene -based  plasmonie  phase  modulator 
seheme  for  THz-band  eommunieation.  The  proposed  modulator  leverages  the  possibility  to  tune 
the  propagation  speed  and,  thus,  output  phase  of  a  SPP  wave  as  its  propagates  over  a  graphene 
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(a)  Constellation  (b)  SER 


Fig.  4:  Performance  of  a  binary  plasmonic  phase  modulator  with  £77^0  =  0.28  eV,^^  !  =  0.13  eV,  (/c=4THz, 
L=2  pm,  t/=90nm). 

layer.  We  have  developed  an  analytieal  model  for  the  modulator  by  starting  from  the  surfaee 
eonduetivity  model  of  graphene  and  the  dispersion  equation  of  SPP  waves  in  gated  graphene 
struetures.  COMSOL  Multi-physies  has  been  utilized  to  validate  the  proposed  model.  We  have  then 
analyzed  the  performanee  of  plasmonie  modulator  in  terms  of  symbol  error  rate  for  the  speeifie  case 
of  a  digital  binary  phase  shift  keying  modulation.  The  results  have  shown  that,  despite  generating  a 
non-uniform  signal  space  constellation,  the  modulated  symbols  are  sufficiently  apart  to  be  easily 
distinguishable.  This  highlights  the  potential  of  the  proposed  approach  to  enable  practical  wireless 
communication  systems  in  the  THz  band. 
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4  GRAPHENE-BASED  PLASMONIC  NANO-ANTENNA  ARRAYS 


4.1  Motivation 

After  modulation,  the  propagating  SPP  wave  needs  to  be  radiated.  As  first  shown  in  [53]  by  using 
the  transmission  line  formalism,  a  graphene -based  antenna  just  one  mierometer  long  and  several 
hundreds  of  nanometers  wide  can  be  designed  to  resonate  in  the  THz-band.  This  frequency  is  two 
orders  of  magnitude  lower  than  the  resonance  frequency  of  a  metallic  antenna  with  comparable 
dimensions.  Nevertheless,  there  are  several  limitations  in  single- antenna  systems  that  motivate  the 
development  of  graphene -based  nano-antenna  arrays  (Fig.  5)  [25].  An  array  configuration  would 
simultaneously  solve  some  of  the  inherent  roadblocks  to  THz  communication  while  opening  up 
many  new  possibilities.  On  the  one  hand,  despite  their  individual  efficiency,  single  nano-antennas 
suffer  from  limited  available  output  power  due  to  their  very  small  size  (a  few  muw?).  On  the  other 
hand,  the  THz  band  itself  introduces  significant  losses  in  the  form  of  spreading  losses  and  attenuation 
due  to  molecular  absorption.  A  THz  band  array  can  alleviate  these  problems  by  producing  huge 
increases  in  gain  due  to  power  amplification  and/or  beamforming.  In  addition,  large  arrays  introduce 
the  possibility  to  take  advantage  of  Multiple-Input  and  Multiple- Output  (MIMO)  communication 
schemes  for  both  beamforming  and  spatial  multiplexing. 


Fig.  5:  Graphene-based  plasmonic  nano-antenna  array. 

In  the  design  of  a  plasmonic  nano-antenna  array,  there  are  two  equally  important  challenges. 
First,  the  individual  nano-antenna  that  serves  as  the  basic  element  of  the  array  must  be  designed. 
This  has  already  been  the  subject  of  multiple  works,  starting  with  [53].  Since  then,  there  have 
been  many  studies  concerning  graphene-based  nano-antennas  which  investigate  variations  such  as 
electrically  tunable  [49],  reconfigurable  [50],  and  beamforming  antennas  [54].  The  second  challenge 
is  to  determine  how  these  individual  elements  should  be  placed  relative  to  each  other.  Here  mutual 
coupling  concerns  become  just  as  important  as  individual  design,  as  this  determines  how  closely 
antennas  can  be  placed  and  subsequently  determines  the  size  of  the  total  array.  Some  works  have 
been  done  in  this  area,  but  they  are  not  sufficient  or  applicable  to  this  case.  For  example,  in  [55] 
coupling  is  considered  between  graphene  nanoribbons,  but  it  does  not  consider  how  the  distance 
between  elements  affects  coupling  strength  and  how  it  impacts  system  performance.  Analysis  is 
needed  which  considers  coupling  among  nano-antennas. 

As  part  of  this  project,  we  have  developed  the  fundamental  theory  of  graphene -based  plasmonic 
nano-antenna  arrays  for  THz  communications.  In  particular,  we  have  first  analyzed  the  mutual 
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coupling  between  graphene-based  nano-antennas  (Sec.  4.2.1).  More  specifically,  we  have  developed 
a  mathematical  analysis  along  with  numerical  simulations  of  a  two-element  array  in  order  to 
understand  how  the  strength  of  mutual  coupling  changes  with  respect  to  the  separation  between 
plasmonic  elements.  The  second  contribution  has  been  an  investigation  of  the  performance  of 
nano-antenna  arrays  in  terms  of  the  achievable  gain  and  directivity  (Sec.  4.3.1).  In  addition  to  the 
ideal  gain,  the  impact  of  mutual  coupling  on  the  gain  has  been  considered  in  our  analysis.  As  before, 
the  developed  models  have  been  validated  by  means  of  electromagnetic  simulations  with  COMSOL 
Multi-physics  (Sec.  4.4). 


4.2  Methods,  Assumptions  and  Procedures 

4.2.1  Mutual  Coupling  Between  Plasmonic  Nano-Antennas 


Coupled  Mode  Theory  The  coupling  between  two  resonant  modes,  such  as  those  excited  in 
resonant  plasmonic  nano-antennas,  can  be  described  mathematically  using  coupled  mode  theory. 
A  full  coupled  mode  model  of  the  nano-antenna  consists  of  a  resonator  with  terms  for  conduction 
losses,  radiation  losses,  and  incident  power  from  outside  waves.  From  [55],  the  amplitude  of  the 
fields  in  two  such  nano-antennas  is  given  by: 


dd\ 

dt 

dd2 

dt 


=  {icoi  -  n  -  ri)di  -I-  ikd2  +  +  iyfne  (21) 

=  (ia)2  -  72  -  r2)d2  -I-  ikdi  +  i^[r2{s+  +  (22) 


where  a\  and  a2  refer  to  the  amplitudes  and  oji  and  a>2  are  the  natural  frequencies  [56].  In  addition, 
7  is  the  radiative  loss,  T  is  the  conduction  loss,  k  is  the  coupling  coefficient,  s+  is  an  incoming  plane 
wave,  and  di,2  =  While  this  is  physically  the  case,  when  the  goal  is  to  investigate  the  near 

field  coupling  the  model  is  greatly  simplified  by  considering  only  the  case  of  two  lossless  resonators. 
This  is  possible  because  the  radiative  and  conductive  losses  are  purely  real,  and  therefore  do  not 
affect  the  resonance  conditions  induced  by  the  reactive  fields.  In  addition,  when  nano-antenna 
separation  distances  are  small  enough  the  near  field  coupling  effects  dominate  and  the  coupling 
due  to  re -radiation  of  energy  can  be  assumed  to  be  zero  [55].  Therefore,  we  assume  that  T  =  0, 
7  =  0,  and  5+  =  0,  and  that  the  resonators  are  excited  by  a  time  varying  function.  In  this  case,  the 
equations  for  the  resonators  are  given  as: 


da\ 

—  =  ]a)\a\  +  kna2, 
dt 

da2 

—  =  jco2a2  +  k2\ai. 
dt 


(23) 

(24) 


The  theory  is  valid  as  long  as  the  coupling  perturbations  can  be  considered  linear,  which  is  true  if 
the  coupling  terms  k\2a2  and  k2\ai  are  much  smaller  than  ja)\a\  and  ja)2a2.  Due  to  the  symmetry 
of  the  problem,  as  well  as  energy  conservation  constraints  discussed  in  [56],  kn  must  be  equal  to 
k2i .  By  taking  this  system  of  linear  equations,  one  can  solve  for  the  eigenfrequencies,  resulting  in: 


to  = 


tO\  +  002 


toi  -  a»2\2 


+  |ki2p. 
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(25) 


where  all  the  terms  have  been  previously  defined. 

Physically  the  eigenfrequencies  can  be  understood  as  the  natural  resonant  frequencies  of 
the  two  nano-antennas.  Because  a  nano-antenna  is  a  passive  component,  if  the  natural  resonance 
frequency  changes  the  result  is  not  that  the  antenna  will  change  operating  frequency,  but  rather  that 
it  will  now  only  radiate  efficiently  at  the  new  resonances.  In  a  system  where  the  nano-antennas  are 
already  being  exciting  by  a  driving  frequency,  a  change  in  the  natural  resonance  will  manifest  itself 
as  a  change  in  the  input  impedance  of  the  system.  A  typical  array  will  be  excited  with  the  same 
frequency  in  each  antenna.  In  this  case,  oJi  =  0)2  =  ojq,  and  the  equation  for  the  eigenfrequencies 
simplifies  to 

oj  =  ojQ  ±  k,  (26) 

where  k  is  the  coupling  coefficient. 

4.3  Coupling  Coefficient 

The  coupling  coefficient  is  the  result  of  the  non-radiating  near  fields  in  the  antenna  and  in  conventional 
microstrip  antenna  theory  it  is  derived  from  Maxwell’s  equations.  There  have  been  multiple  attempts 
to  model  the  mutual  coupling  between  elements  by  using  the  transmission  line  model  [57],  the 
cavity  model  [58]  [59]  [60],  or  the  method  of  moments  for  microstrip  antennas.  While  these  studies 
have  had  success  in  predicting  coupling  in  metallic  microstrip  antennas,  they  have  drawbacks  which 
limit  their  use  in  studying  coupling  between  plasmonic  nano-antenna  arrays.  The  problem  is  that 
most  studies  neglect  surface  waves  to  simplify  the  model.  In  a  plasmonic  nano-antenna,  surface 
waves  are  critical  to  the  operation  of  the  antenna,  and  in  fact  are  the  only  kind  of  waves  that  can  be 
supported  [19].  Even  when  studies  account  for  surface  waves,  these  surface  waves  are  fundamentally 
different  from  SPP  waves.  The  surface  waves  on  metallic  antennas  occur  at  dielectric -dielectric 
boundaries  and  propagate  over  the  entire  dielectric  surface.  These  waves  detract  from  the  total  power 
available  to  be  radiated.  In  contrast,  a  plasmonic  nano-antenna  relies  on  SPP  waves  for  generating 
the  electric  field  and  are  only  supported  on  the  boundary  between  a  dielectric  and  a  metal,  which 
means  they  can  only  exist  where  the  patch  is  present.  Because  of  these  differences,  the  models  used 
for  predicting  the  coupling  coefficient  in  metallic  antenna  arrays  cannot  be  assumed  to  work  for 
plasmonic  arrays.  Further  analysis  is  needed  which  takes  into  account  the  properties  of  SPP  waves. 

Because  the  wavelength  in  a  plasmonic  nano-antenna  is  confined  to  a  much  smaller  length 
than  the  radiated  free  space  wavelength,  it  is  expected  that  the  mutual  coupling  will  only  start  to 
become  a  major  factor  at  distances  related  to  the  plasmonic  wavelength,  and  consequently  allow 
array  elements  to  be  separated  much  less  than  would  be  required  in  a  perfect  electric  conducting 
antenna.  Taking  this  into  consideration,  we  predict  that  the  coupling  coefficient  can  be  approximated 
by  an  exponential  function  given  as 

k  =  a(Ooe~‘^^,  (27) 

where  ojq  is  the  resonant  frequency  without  coupling,  a  and  j3  are  tuning  constants,  and  d  is  the 
separation  distance  between  elements. 

4.3.1  Performance  of  Plasmonic  Nano-antenna  Arrays 

In  the  design  of  an  array,  it  is  important  to  consider  the  overall  system  performance  and  how  it 
is  impacted  by  coupling  and  other  factors.  An  array  constructed  using  graphene  nano-antennas 
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will  differ  in  two  important  aspects  from  a  conventional  array.  First,  as  previously  discussed,  the 
coupling  is  expected  to  be  based  mainly  on  the  physical  size  of  the  nano-antenna.  For  a  plasmonic 
nano-antenna  where  there  is  a  high  conhnement  factor,  near  held  coupling  will  only  be  an  issue  at 
distances  much  less  than  the  free  space  wavelength.  On  one  hand,  this  will  allow  for  array  elements 
to  be  placed  much  closer  in  proximity  to  each  other,  resulting  in  a  high  density  of  elements  in  the 
available  area.  On  the  other  hand,  this  limits  the  beamforming  abilities,  which  require  distance 
comparable  to  the  free  space  wavelength  in  order  to  achieve  constructive  superposition  of  held 
amplitudes.  Second,  the  possibility  to  use  one  plasmonic  signal  source  per  antenna  to  power  every 
element  individually  means  that  an  increase  in  output  power,  and  consequently  gain,  can  be  achieved 
independently  of  phase  variation  or  beamforming.  These  differences  are  illustrated  by  considering 
the  array  factor  of  a  uniform  square  planar  array  with  N  powered  elements  per  side,  separated 
by  a  distance  of  A,  which  is  given  in  [61]  by 


^sin 

(wt 

)' 

(  sin 

)' 

^  sin 

(fe) 
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/ 

4^x  — 


2n 

kd  sin  6  cos  0  -  —  sin  6q  cos 
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iky  = 


kd  sin  0  sin  0 


In 
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sin  Oq  sin  (pQ, 


(28) 

(29) 

(30) 


where  6  refers  to  the  elevation  angle,  (p  is  the  azimuth,  (0o,  is  the  beam  pointing  direction,  k  =  ^ 
is  the  wavenumber,  and  d  is  the  uniform  distance  between  elements.  In  a  conventional  array 
is  normalized  with  respect  to  N  to  account  for  the  total  power  being  split  between  N  antennas, 
and  to  better  illustrate  radiation  pattern  changes  from  held  multiplication.  However,  here  it  is  left 
unnormalized  with  respect  to  N  due  to  the  fact  that  the  nano-antennas  will  be  powered  independently 
and  contribute  to  the  overall  gain.  In  addition,  for  a  plasmonic  nano-antenna  array  the  distance 
between  the  elements  will  be  the  plasmonic  wavelength  rather  than  some  multiple  of  the  freespace 
wavelength.  In  this  case,  we  replace  d  with  the  plasmonic  wavelength  Tspp,  and  the  array  factor 
becomes 


^T{e,  (p) 

ipx 


^sin 

)' 

^sin  {n'^^ 

)' 

^  sin 

(t) 

/ 

2n  2n 

—  sin  6  cos  (p - sin  6q  cos  (po, 

y  y 


2n  2n 

—  sin  y  sin  0 - sin  6q  sin  cpQ, 

y  y 


(31) 

(32) 

(33) 


where  y  is  the  conhnement  factor  and  is  equal  to  ^ .  When  using  this  formula  it  is  clear  that 
the  higher  the  conhnement  factor  becomes,  the  lower  the  directivity  will  be,  although  it  still  has 
a  minimum  value  of  1 .  However,  analysis  using  this  formula  suggests  that  even  with  decreased 
directivity  due  to  a  high  conhnement  factor,  or  equivalently,  close  element  spacing,  the  attainable 
gain  is  still  high  due  to  power  multiplication. 
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Table  1 :  Curve  fitted  Parameters  for  k 


Frequency 

1  THz 

5  THz 

8  THz 

a 

0.02563 

0.02158 

0.02188 

P 

lel0““ 

lelO-^^ 

lel0-“ 

4.4  Results  and  Discussions 

In  this  section,  the  analytical  model  for  the  mutual  coupling  coefficient  is  validated,  and  the  nano¬ 
antenna  array  performance  claims  are  verified  by  means  of  full-wave  electromagnetic  simulations.  In 
order  to  numerically  find  the  coupling  coefficient,  a  two  element  array  is  simulated  using  COMSOL 
Multiphysics.  Graphene  is  modeled  as  a  transition  boundary  condition  with  complex-valued  dynamic 
conductivity  given  by  the  model  in  [27],  where  Tg  =  0.5  ps  is  the  relaxation  time  of  electrons  in 
graphene,  T  =  300  K  is  the  temperature,  and  with  Ef  values  from  0.1  to  1.25  eV,  where  Ef  is  the 
Fermi  energy  of  the  graphene  patch.  These  values  are  based  on  analysis  of  Raman  spectra  obtained 
from  CVD-grown  graphene.  The  graphene  layer  rests  on  top  of  a  90-nm-thick  Si02  dielectric 
{er  =  4),  which  separates  the  graphene  from  the  ground  plane.  This  thickness  is  chosen  because 
it  maximizes  visual  detection  of  graphene  on  Si02.  The  antenna  is  fed  with  a  lumped  port  that 
connects  the  graphene  layer  to  the  ground  plane  on  one  side.  A  perfectly  matched  layer  and  scattering 
boundary  are  used  to  accurately  approximate  an  infinite  space.  The  antennas  are  meshed  with  a 
resolution  of  XsppjS.  To  simulate  the  array  performance  measures  of  gain  and  directivity,  the  same 
parameters  are  used  to  construct  a  uniform  planar  three  by  three  array  for  a  total  of  nine  elements. 

4.4.1  Simulation  Validation  of  Mutual  Coupling 

The  effect  of  mutual  coupling  on  resonant  frequency  was  found  through  the  measurement  of  the 
input  impedance  at  the  ports.  As  previously  discussed,  the  frequency  at  which  the  imaginary  part  of 
the  input  impedance  disappears  is  considered  to  be  the  resonant  frequency.  The  mutual  coupling 
dependency  on  the  separation  of  the  elements  was  found  by  performing  a  frequency  sweep  for 
different  separation  distances  for  three  different  frequencies  that  are  representative  of  the  THz 
band.  The  coupling  coefficient  (27)  was  also  plotted  using  curve  fitting  tools  to  match  the  tuning 
parameters.  As  shown  in  Fig.  6,  the  change  in  frequency  as  a  function  of  separation  distance  is 
comparable  to  the  plasmonic  wavelength,  which  in  this  case  is  2  pm  for  both  simulations.  It  also 
shows  that  the  matched  coupling  coefficient  equations  are  well  matched  to  their  respective  coupling 
sweeps.  Table  1  summarizes  the  matched  coefficients  for  each  frequency. 

4.4.2  Numerical  Analysis  of  Array  Performance 

In  addition  to  studying  the  change  in  frequency  in  a  two  element  array,  the  overall  array  performance 
and  the  impact  of  mutual  coupling  on  this  performance  was  simulated.  The  simulation  used  a 
nine  element  array,  and  directivity  and  gain  are  calculated  from  radiated  power  values  computed 
by  COMSOL.  The  directivity  as  plotted  in  Fig.  7a  takes  into  account  the  combined  effects  of  the 
array  directivity  and  the  single  nano-antenna’s  directivity,  while  neglecting  power  amplification. 
What  can  be  seen  from  the  graph  is  that  the  directivity  when  the  antenna  separation  distance  is 
comparable  to  Tspp  is  predominately  due  to  the  directivity  of  the  single  nano-antenna,  which  is  about 
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(a)  Plasmonic  Array  Resonant  Frequency  at  1  THz 


x10'2 


(b)  Plasmonic  Array  Resonant  Frequency  at  5  THz 


(c)  Plasmonic  Array  Resonant  Frequency  at  8  THz 


Fig.  6:  Resonance  frequency  as  a  function  of  the  separation  distance  between  resonant  elements  for  different 
antenna  resonant  frequencies. 

6.5  dB.  The  added  direetivity  from  the  array  due  to  eonstruetive  interferenee  does  not  oeeur  until 
the  separation  approaehes  Ajl.  In  Fig.  7b,  the  power  multiplieation  gain  of  the  array  over  that  of  a 
single  nano-antenna  is  plotted.  This  gain  in  power  is  shown  to  be  equal  in  magnitude  to  the  square 
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of  the  array  factor,  which  for  an  array  of  nine  elements  is  equivalent  to  81,  or  about  19  dB.  In  both 
figures  it  is  clear  that  mutual  coupling  effects  cause  a  large  drop  in  both  directivity  and  power  gain, 
but  only  for  very  close  separation  distances. 


(a)  Plasmonic  Array  Directivity  at  8  THz 


(b)  Plasmonic  Array  Power  Gain  at  8  THz 

Fig.  7:  Directivity  and  power  gain  of  a  nine  antenna  uniform  planar  array  as  functions  of  the  separation 
between  elements. 

Finally,  in  Fig.  8,  the  total  gain  of  a  plasmonic  array  is  calculated  based  on  the  number  of 
antennas  in  the  array  for  different  separation  distances.  The  total  gain  shown  here  represents  the 
maximum  radiated  power  of  a  plasmonic  array  compared  to  that  of  a  single  isotropic  nano-antenna 
and  includes  the  directivity  of  a  graphene  nano-antenna,  the  directivity  of  the  array,  and  the  power 
multiplication  from  multiple  powered  antennas.  It  is  important  to  note  that  even  the  lowest  separation 
distance  shown  here  is  large  enough  that  mutual  coupling  effects  are  not  a  concern. 

4.4.3  Discussion 

As  part  of  this  effort,  we  have  analyzed  the  performance  of  graphene-based  plasmonic  nano-antenna 
arrays  at  THz  frequencies,  by  taking  the  impact  of  mutual  coupling  into  account.  Coupled  mode 
theory  has  been  utilized  to  model  the  mutual  coupling  between  two  neighboring  plasmonic  nano¬ 
antennas  and  the  coupling  coefficient  has  been  defined.  The  resulting  gain  of  an  active  plasmonic 


Approved  for  Public  Release;  Distribution  Unlimited 

20 


400  600  800 

Number  of  Elements 


1200 


Fig.  8:  Total  gain  of  a  graphene-based  plasmonic  nano-antenna  array,  as  a  function  of  the  number  of  elements. 


nano-antenna  array  has  been  derived,  and  finite-element  simulations  have  been  utilized  to  both 
validate  the  models  and  numerieally  explore  the  performanee  of  plasmonie  nano-antenna  arrays.  It 
has  been  shown  that  near  field  eoupling  between  antennas  beeomes  negligible  at  separation  distanees 
mueh  less  than  the  free  spaee  wavelength.  Consequently,  a  very  large  number  of  elements  ean  be 
integrated  in  plasmonie  nano-antenna  arrays  with  a  resulting  very  small  footprint  and,  potentially, 
lead  to  practical  THz  communication  systems. 
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5  MULTI-USER  INTERFERENCE  MODELING  IN  THE  THZ  BAND 

5.1  Motivation 

As  mentioned  in  the  introduction,  the  THz  band  provides  wireless  communication  devices  with  an 
unprecedentedly  large  bandwidth,  ranging  from  several  tens  of  GHz  up  to  a  few  THz  [23, 62].  The 
main  phenomenon  affecting  the  propagation  of  THz-band  signals  is  the  absorption  by  water  vapor 
molecules.  For  communication  distances  below  1  m,  where  the  number  of  molecules  found  along 
the  path  is  small,  the  THz  band  behaves  as  a  single  transmission  window  several  THz  wide.  In  this 
scenario,  the  transmission  of  one-hundred-femtosecond-long  pulses  by  following  an  on-off  keying 
modulation  spread  in  time  has  been  recently  proposed  [63].  The  power  spectral  density  (PSD) 
of  such  pulses,  which  can  be  generated  and  detected  by  means  of  different  electronic,  photonic 
and  plasmonic  devices,  has  its  main  frequency  components  under  4  THz.  By  utilizing  an  on-off 
keying  modulation,  the  transmitter  can  reduce  its  energy  consumption  by  remaining  silent  during 
the  transmission  of  logical  “Os",  and  by  spreading  the  transmission  of  the  pulses  in  time,  concurrent 
transmissions  can  be  multiplexed  in  time. 

The  transmission  of  ultra-short  pulses  minimizes  the  probability  of  collisions  due  to  the 
very  small  time  that  the  channel  is  occupied  by  each  user.  However,  given  that  many  of  the 
envisioned  applications  of  THz-band  communications  involve  very  large  node  densities  (e.g., 
wireless  nanosensor  networks  or  networks  on  chip),  multi-user  interference  results  unavoidably. 
Many  stochastic  models  of  interference  have  been  developed  to  date  [64-66].  However,  these  models 
do  not  capture  the  peculiarities  of  the  THz  band  channel,  such  as  the  molecular  absorption  loss. 
More  recently,  new  interference  models  for  THz-band  communication  have  been  developed  [63, 67]. 
However,  these  models  do  not  take  into  account  the  specific  waveforms  that  are  being  transmitted 
(e.g.,  Gaussian-shaped  pulses),  and,  in  addition,  no  experimental  validation  for  any  of  the  models 
exists  to  date. 

As  part  of  this  effort,  we  have  developed  a  stochastic  model  of  multi-user  interference  in 
pulse-based  THz-band  communications  and  experimentally  validated  it.  This  model  is  developed  by 
considering  the  fact  that  the  interference  power  at  the  receiver  is  not  a  combination  of  the  received 
powers  from  the  individual  nodes,  rather  the  power  of  the  combination  of  the  signal  amplitudes. 
For  this,  first,  we  obtained  the  probability  density  function  (PDF)  of  the  interference  generated  by 
one  interfering  node  at  the  receiver,  starting  from  the  PDFs  of  the  pulse  received  energy  and  the 
PDF  of  the  pulse  shape.  Then,  we  extended  this  model  to  account  for  multiple  nodes  which  can 
constructively  or  destructively  interfere.  Finally,  we  validated  the  model  by  comparing  the  results  to 
experimental  measurements.  Next,  in  Sec.  5.2.1,  we  describe  the  system  model.  The  stochastic 
model  of  multi-user  interference  is  developed  in  Sec.  5.2.2.  Finally,  we  describe  our  experimental 
setup,  validate  the  analytical  model  and  provide  numerical  results  in  Sec.  5.3. 

5.2  Methods,  Assumptions  and  Procedures 
5.2.1  System  Model 

To  model  the  multi-user  interference,  we  first  need  to  take  into  account  both  the  spatial  distribution 
of  the  nodes  as  well  as  their  temporal  activity.  In  this  section,  we  describe  the  system  model  utilized 
in  the  paper. 
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Network  Topology  We  consider  that  nodes  are  randomly  distributed  in  space  by  following  a 
spatial  Poisson  process  with  rate  A.  Without  losing  generality,  we  consider  that  the  intended  receiver 
is  located  at  the  center  of  a  disc  of  radius  a  and  area  A  =  na^.  The  probability  of  finding  k  nodes  in 
the  disc  is  given  by 

(ATra^)^  .2 

Pike  A)  =  ^  .  (34) 

k\ 

Under  the  Poisson  assumption,  the  locations  of  the  nodes  follow  independent  and  identically 
distributed  uniform  distributions.  If  d  is  the  distance  to  the  origin  from  a  point  that  is  uniformly 
distributed  in  A,  then  the  PDF  of  the  distance  distribution  D  is  given  by 


(2d 


fDid)  = 


for  0  <  d  <  a 
otherwise. 


(35) 


Physical  and  Link  Layers  We  consider  that  nodes  communicate  by  utilizing  TS-OOK  [63].  In 
this  scheme,  a  symbol  “1"  is  sent  by  transmitting  a  hundred-femtosecond  long  pulse  and  a  symbol 
“0"  is  sent  by  silence,  i.e.,  a  node  does  not  transmit  anything  (Fig.  9).  The  time  between  symbols 
is  much  longer  than  the  symbol  duration  Tp,  i.e.,  /3  =  T^/Tp  »  1. 


x10® 


Time[s]  xIO'^^ 


Fig.  9:  Illustration  of  a  TS-OOK  transmission. 


The  pulses  in  TS-OOK  are  usually  modeled  as  Gaussian-shaped.  More  specifically,  the 
radiated  pulses  s  (t)  resemble  the  first  time  derivative  of  a  Gaussian  pulse. 


sit)  = 


dx  it) 


dt 


- 2 — -^(0^  wrth 


cr^ 


X  it)  =  ao 


_(tztr 
e  2o-2 


V2^i 


(36) 

(37) 


no¬ 


where  ao  is  a  normalizing  constant  to  adjust  the  pulse  total  energy,  cr  is  the  standard  deviation  of 
the  Gaussian  pulse  in  seconds,  and  yu  is  the  center  of  the  pulse  in  time,  given  in  seconds.  The  total 
length  of  the  radiated  pulses  ^  it)  is  approximately  Tp  =  lcr. 

TS-OOK  enables  robust  and  concurrent  communication  among  nodes.  In  the  envisioned 
scenarios,  nodes  can  start  transmitting  at  any  time  without  being  synchronized  or  controlled  by  any 


Approved  for  Public  Release;  Distribution  Unlimited 

23 


type  of  network  eentral  entity.  However,  due  to  the  faet  that  the  time  between  transmissions  Ts  is 
mueh  longer  than  the  pulse  duration  Tp,  several  nodes  ean  eoneurrently  use  the  ehannel  without 
neeessarily  affeeting  eaeh  other.  In  addition,  the  very  short  symbol  duration  Tp  makes  eollisions 
between  symbols  highly  unlikely.  However,  as  the  number  of  eommunieating  nodes  inereases, 
interferenee  beeomes  a  problem,  as  we  show  in  the  next  seetion. 

5.2.2  Interference  Model 

In  the  existing  interferenee  models,  the  total  interferenee  power  at  the  reeeiver  is  obtained  as 
the  addition  of  the  interferenee  power  from  eaeh  node.  However,  by  eonsidering  the  reeeived 
signal  power  instead  of  the  reeeived  signal  amplitude,  the  faet  that  interferenee  ean  be  either 
eonstruetive  or  destruetive  is  effeetively  negleeted.  This  usually  results  in  unrealistieally  large  values 
of  interferenee,  speeially  when  the  transmitted  signals  have  large  fluetuations  between  positive  and 
negative  amplitude  values.  Our  analytie  model  overeomes  this  limitation  by  taking  into  aeeount  both 
the  signal  waveform  as  well  as  the  THz-band  ehannel.  We  taylor  our  model  to  Gaussian-shaped 
pulses,  but  the  same  methodology  eould  be  utilized  for  any  other  waveform. 

In  order  to  eharaeterize  the  interferenee,  we  observe  the  signals  at  the  reeeiver.  For  the  system 
model  deseribed  in  See.  5.2.1,  the  signal  y  (t)  at  the  reeeiver  ean  be  written  as 

u 

y{t)  =  ^  (5  (i  -  Tu))  *  h  (t,  du)  +  n(t),  (38) 

U=1 

where  U  is  the  total  number  of  nodes  in  the  network,  s  (t  -  Tu)  is  the  transmitted  signal  along  with 
the  delay  r^,  h  (t,  du)  is  the  system  impulse  response,  *  denotes  eonvolution,  and  n  (t)  is  the  noise  at 
the  reeeiver.  The  system  impulse  response  h  {t,  du)  eaptures  the  impaet  of  the  THz-band  ehannel, 
and  is  obtained  as  in  [63]. 

If  we  eonsider  the  signal  from  user  m  =  1  as  the  intended  signal,  the  total  interferenee  is  given 
by 

u 

i  (0  =  ^  iu  it) ,  (39) 

u=2 

where  (t)  is  the  interferenee  at  reeeiver  due  to  eaeh  node  u.  By  foeusing  on  one  speeihe  symbol, 
without  losing  generality,  iu  (t)  ean  be  written  as 

iu  (t)  =  ^/^P  {t  -  Tf,„) ,  (40) 

where  Cp^u  stands  for  the  energy  of  the  reeeived  pulse,  p  (t)  is  the  shape  of  the  pulse  with  unitary 
energy,  and  Tt^u  refers  to  the  total  delay  (initial  delay  and  propagation  delay). 

Given  that  nodes  are  randomly  distributed  in  spaee  and  uneoordinatedly  transmitting  in  time, 
the  PDF  of  the  interferenee  distribution  /„  generated  by  a  node  u  ean  be  modeled  as  a  funetion  of 
two  random  variables,  namely,  the  energy  of  the  pulse,  with  distribution  Ep,  and  the  shape  of  the 
pulse,  with  distribution  P. 

On  the  one  hand,  the  energy  of  the  pulse  depends  on  the  transmitted  pulse  energy,  whieh  we 
eonsider  a  parameter  in  our  analysis,  and  the  distanee  between  the  nodes,  whieh  is  a  random  variable 
D  with  PDF  given  by  (35).  Therefore,  the  PDF  of  Ep  ean  be  derived  from  the  PDF  of  D  as  follows: 
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Fig.  10:  Sinusoidal  approximation  of  the  transmitted  pulse  in  TS-OOK. 


fEp  = 


fDir\e,)) 


(41) 


where  denotes  the  inverse  and  f  denotes  the  first  derivative  of  Cp.  The  PDF  of  Ep  can  only  be 
analytically  obtained  if  Cp  (d)  is  invertible.  For  the  distances  considered  in  our  analysis,  it  can  be 
shown  that  the  received  energy  can  be  approximated  by  the  polynomial 


e 


p 


(42) 


where  ^  and  rj  are  two  constants  which  depend  on  the  specific  channel  molecular  composition  as 
well  as  on  the  energy  and  the  shape  of  the  transmitted  signal.  By  computing  the  inverse  function 
of  (42),  and  combining  it  with  (35)  in  (41),  we  obtain 


fEp  [ep] 
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2+1] 
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for  ^  <  Cp  <  oo 
otherwise. 


(43) 


The  step  by  step  derivation  is  methodological  and  has  been  omitted  due  to  space  constraints.  Finally, 
the  interference  from  a  node  depends  on  the  amplitude  of  the  received  pulse,  not  its  energy.  If  we 
dehne  Cs  =  ^Je^,  the  PDF  of  Es  is  now  given  by 


fEs  {es) 


A^^l^\es\  ~ 

a^T] 
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for  <  es  <  oo 

otherwise. 


(44) 


On  the  other  hand,  in  order  to  capture  the  constructive  or  destructive  addition  of  the  interference 
at  the  receiver,  we  need  to  take  into  account  the  shape  of  the  transmitted  pulse,  which  is  deterministic, 
and  the  time  instant  at  which  it  reaches  the  receiver,  which  we  model  as  a  random  variable  T  with 
uniform  distribution  over  the  duration  between  two  symbols  Ts, 
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(45) 
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As  before,  to  obtain  the  PDF  of  P,  we  need  to  first  obtain  the  inverse  function  of  p  (t),  which  in 
this  case  is  given  by  (37).  However,  this  cannot  be  obtained  in  closed  form  expression.  To  overcome 
this  limitation,  we  approximate  the  transmitted  pulse  with  a  shifted  sinusoidal  pulse  of  same  duration 
Tp,  as  shown  in  Fig.  10.  This  pulse  can  be  expressed  mathematically  as  follows 

=  (46) 


where  r  =  is  the  amplitude  of  the  sinusoid  normalized  to  have  unit  energy.  The  inverse  function 
of  this  new  pulse  can  be  written  as 


Now  the  PDF  of  P  can  be  derived  from  the  PDF  of  T  similarly  as  iox  Ep,  yielding  to 


P  .4.-1  (P 


(47) 


fp{p)  = 


T.nrJl 
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for  -r  <  ;?  <  r 


otherwise. 


(48) 


T 

Note  that  the  above  PDF  integrates  to  #,  not  1.  This  is  because  the  pulse  p  has  nonzero  values  only 

^  s 

T 

within  a  duration  of  Tp.  For  rest  of  the  time  1  -  the  signal  value  is  zero.  Hence,  we  express  the 
PDF  of  P  as  a  mixed  distribution 

fp(p)  = - ^ - -  +  (l-^)<)(rt.  (49) 

where  6  stands  for  the  Dirac-delta  function. 

Since  the  interference  generated  by  one  user  is  the  product  of  two  random  variables  Es  and  P, 
its  PDF  can  be  determined  from  their  PDFs  as  follows: 

flu  (iu)  =  f  fp(p)fEs  (-)  -dp.  (50) 

J -CO  \p  I  P 

By  replacing  the  PDFs  from  (44)  and  (49)  in  (50),  the  PDF  of  lu  is  obtained  (51),  where  F  () ,  B  () , 
and  2^1  0  are  the  standard  Gamma,  Beta  and  Hyper-geometric  functions: 
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(51) 


Since  the  signals  from  individual  nodes  add  up  at  the  receiver  and  are  independent,  the  PDF  of 
the  total  interference  can  be  determined  by  the  convolution  of  the  individual  PDFs: 

fi  (0  =  fh  ih)  *  fh  (h)  *  fh  ih)  flu  iw)  ■  (52) 

To  calculate  the  n-fold  convolution  on  the  right  hand  side,  we  go  into  transform  domain.  We  define 
the  characteristic  function  of  the  interference  from  node  u  as 

=  (53) 

Since  the  interferences  from  nodes  are  independent  and  identically  distributed  random  variables, 
they  have  the  same  characteristic  function.  The  convolution  becomes  a  product  in  the  transform 
domain  and,  thus,  the  characteristic  function  of  the  interference  can  be  written  as 

(p{io)  =  {(pu{ioyf  ■  (54) 

Now  the  PDF  of  the  total  interference  can  be  found  by  taking  the  inverse  transform  of  (p  (<^): 

//(0  =  F-i{<^(m)}.  (55) 

Unfortunately,  this  characteristic  function  cannot  be  analytically  derived  in  closed  form,  and,  thus, 
we  numerically  obtain  it  in  the  next  section. 

5.3  Results  and  Discussions 
5.3.1  Experimental  Setup 

In  order  to  experimentally  emulate  THz  communication  signals  emanating  from  multiple  THz 
nodes,  we  utilize  a  THz  time-domain  spectroscopy  (THz-TDS)  system  combined  with  an  innovative 
interference  generating  setup  (Fig.  11).  In  particular,  a  collimated  3  cm  wide  THz-TDS  beam  is  used 
to  illuminate  a  row  of  metal  rods.  The  metallic  rods  are  1  mm  in  diameter  and  can  effectively  reflect 
the  THz  pulses.  Therefore,  each  rod  effectively  represents  an  interferer  at  a  different  distance. 


Fig.  i  i :  Illustration  of  the  THz  beam  from  the  transmitter  (Tx)  to  receiver  (Rx)  via  reflections  from  the  metal 
rods,  which  emulate  interference  from  multiple  nodes. 


For  the  measurements,  an  apparatus  (Fig.  12)  was  created  to  hold  up  to  20  metal  rods  in  place, 
linearly.  The  THz  transmitter  is  positioned  so  that  the  THz  beam  is  centered  at  position  13,  27.5  cm 
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Fig.  12:  Top  view  of  the  linear  apparatus  that  holds  the  metal  rods  in  place. 

away.  The  receiver  is  set  at  50  degrees  from  the  transmitter  at  a  distance  of  14.5  cm  away  from  a 
metal  rod  in  position  13  (Fig.  11).  Reflected  pulse  obtained  by  the  receiver  shows  the  presence  of 
a  metal  rod.  Data  is  collected  by  taking  scans  of  the  metal  rods  placed  in  specific  positions  and 
obtaining  the  results  in  the  time  domain.  Fig.  13  shows  the  THz  time  domain  waveform  for  the 
reflection  off  of  one  metal  rod  in  position  13. 


Fig.  13:  THz  pulse  from  the  reflection  off  of  one  metal  rod  in  position  13. 


Distance  [m] 


Fig.  14:  Measured  and  approximated  received  pulse  energy  as  functions  of  the  transmission  distance. 

Fig.  14  shows  the  received  energy  versus  distance  curve  for  the  measurements  taken  in 
the  above  described  setup.  A  curve  fitting  was  done  to  determine  the  ^  and  rj  values  from  the 
experimental  data  and  has  been  found  to  be  3.7820e-16  and  0.9357  respectively.  The  value  of  rj 
confirms  the  cylindrical  spreading  of  the  signal  as  it  was  collimated  in  the  experimental  setup.  These 
extracted  parameter  values  have  then  been  used  in  the  analytical  model  to  generate  the  PDF  of  the 
interference  for  different  number  of  nodes  in  the  network.  The  other  parameters  of  the  model  have 
been  chosen  to  match  the  experimental  setup.  For  example,  =  8  ps  has  been  approximated  from 
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the  radiated  pulse  duration,  jS=10  has  been  chosen  to  match  the  observation  period  of  the  measured 
signal.  Also,  a  =  ' 
pulse  unit  energy. 


/  '7 

signal.  Also,  a  =  0.8  m  was  chosen  to  reflect  the  experimental  setup  and  r  =  to  ensure  the 


5.3.2  Single  Interferer  Model  Validation 


Fig.  15:  Interference  distribution  of  a  single  interferer. 


Fig.  15  shows  the  the  probability  distribution  function  of  the  interference  for  both  experimental 
and  analytical  case  when  there  is  a  single  interferer  in  the  network.  As  expected,  most  of  the 
probability  mass  is  concentrated  at  zero  as  shown  by  the  impulse  at  the  center.  The  impulse  was  cut 
to  show  the  graph  of  the  distribution  clearly.  It  can  also  be  seen  from  the  figure  that  the  experimental 
data  has  a  lot  of  samples  near  zero  whereas  the  analytical  model  doesn’t.  The  reason  for  this  is  that 
there  is  inevitable  noise  in  the  experiment  and  we  haven’t  considered  the  noise  in  our  analytical 
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Fig.  16:  Interference  distribution  of  two  interferers. 


model.  Also,  we  couldn’t  take  measurements  below  5  cm  distance  between  the  transmitter  and 
receiver  due  to  the  experimental  set  up  limitations.  Nevertheless,  the  distribution  from  the  analytical 
model  and  experimental  data  follows  similar  trend.  It  can  reasonably  be  said  that,  as  the  number  of 
measurements  is  increased  and  the  distances  are  picked  randomly,  the  experimental  PDF  will  more 
closely  follow  the  analytical  model. 
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5.3.3  Multiple  Interferer  Model  Validation 


Fig.  16  compares  the  interference  PDFs  from  analytical  model  and  experimental  data  for  two 
interferers  in  the  network.  It  can  be  seen  that  they  match  much  better  than  the  single  interferer 
case.  The  reason  is  that  the  measurements  were  taken  by  putting  two  rods  at  random  locations  in 
the  linear  apparatus  which  essentially  ereated  randomness  in  the  distanee  between  transmitter  and 
the  reeeiver.  Fig.  17  shows  that  the  interferenee  spreads  away  from  the  0  value  as  the  number  of 
interferer  inereases  in  the  network.  In  addition,  the  height  of  the  impulse  is  also  deereased,  redueing 
the  ehanee  of  zero  interferenee.  For  large  number  of  interferers,  the  PDF  eonverges  to  a  Gaussian 
PDF  as  expeeted  from  the  Central  Limit  theorem,  and  originally  shown  for  TS-OOK  in  [63]. 


Fig.  17:  Illustration  of  PDF  converging  to  Gaussian  with  increasing  number  of  nodes. 


5.3.4  Discussion 

An  stochastic  model  for  short  range  pulse-based  THz-band  communication  network  has  been 
presented  considering  the  fact  that  interference  occurs  with  the  amplitude  of  the  signal,  not  the 
power.  The  PDF  of  the  interference  has  been  analytically  deduced  from  the  PDFs  of  the  energy 
and  the  pulse.  In  addition,  it  has  been  compared  with  experimental  data  to  validate  the  model. 
The  results  show  similarity  in  spite  of  the  limitation  we  had  in  the  experimental  set  up  and  noise 
modeling.  It  has  also  been  shown  that  the  PDF  converges  to  a  Gaussian  distribution  as  the  number 
of  nodes  increases  in  the  network. 
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6  TERAHERTZ  SYMBOL  DETECTION  AND  PHYSICAL  LAYER  SYNCHRONIZA¬ 
TION 

6.1  Motivation 

Independently  of  interferenee,  one  of  the  main  ehallenges  with  ultra-high- speed  pulse -based 
eommunieations  is  the  need  for  tight  symbol  synehronization  between  transmitter  and  reeeiver. 
Physieal-layer  synehronization  affeets  the  sueeessful  symbol  deteetion  probability  and  impaets 
the  overall  link  performanee.  Unfortunately,  we  eannot  direetly  reuse  the  existing  solutions  for 
similar  pulse-based  eommunieation  systems,  sueh  as  Impulse-Radio  Ultra  Wide  Band  (IR-UWB) 
eommunieations  or  Free  Spaee  Optieal  systems  (FSO)  [68-70].  One  of  the  main  reasons  for  this  is 
the  very  short  duration  of  THz  pulses,  whieh  require  the  use  of  very  high-speed  Analog-to-Digital 
Converters  (ADCs)  for  synehronization  and  symbol  deteetion.  The  fastest  existing  ADC  to  date  ean 
only  sample  at  rates  below  100  Giga-Samples-per-seeond  (GSas)  [71],  mueh  below  the  Nyquist  rate 
for  THz  signals.  Furthermore,  its  size  and  power  eonsumption  make  it  inadequate  for  nano-deviees. 
In  addition  to  the  laek  of  ADCs,  the  loeal  eloek  [72]  at  different  nano-deviees  might  oseillate  at 
slightly  different  frequeneies,  whieh  ean  result  in  a  signifieant  eloek  skew  between  the  transmitter 
and  the  reeeiver. 

As  part  of  this  effort,  we  have  developed  a  new  synehronization  seheme  for  pulse -based 
THz -band  eommunieation,  whieh  is  based  on  dynamieally  time- shifting  the  signal  at  the  reeeiver 
by  means  of  a  voltage-eontrolled  delay  (VCD)  line  [73]  and  adapting  the  observation  window  of  a 
Continuous-Time  Moving- Average  (CTMA)  symbol  deteetor  [74].  The  fundamental  idea  is  that, 
during  an  initial  synehronization  preamble,  a  synehronization  bloek,  whieh  is  eomposed  by  a  beneh 
of  VCD  lines  and  CTMA  deteetors,  is  used  to  iteratively  determine  the  symbol  start  time  and 
observation  window  length  that  is  then  used  for  the  deteetion  of  the  data  symbols  in  the  paeket. 

Next,  we  first  deseribe  and  analyze  the  aeeuraey  of  the  synehronization  bloek  and  provide 
elosed-form  expressions  for  the  number  of  bits  in  the  preamble  neeessary  to  aehieve  synehronization, 
whieh  depends  on  the  number  of  elements  in  the  synehronization  bloek  as  well  as  the  eloek 
skew  between  transmitter  and  reeeiver  (See.  6.2.1).  Then,  we  analyze  the  sueeessful  symbol 
deteetion  probability  of  the  CTMA  symbol  deteetor  as  a  funetion  of  the  observation  window  length 
determined  from  the  synehronization  bloek  (See.  6.2.2).  In  addition,  we  investigate  the  impaet  of 
the  synehronization  and  symbol  deteetion  bloeks  on  the  aehievable  link  throughput,  and  identify  the 
optimal  parameter  values  for  their  design  (See.  6.2.3).  Moreover,  we  validate  and  test  the  seheme 
with  experimental  data  obtained  with  a  THz  time-domain  speetroseopy  platform,  and  provide 
extensive  numerieal  results  to  support  our  developed  solution  (See.  6.3). 

6.2  Methods,  Assumptions  and  Procedures 

6.2.1  System  Model  and  Synchronization  Algorithm 

In  our  scenario,  we  consider  that  nodes  communicate  in  an  uncoordinated  fashion  and  by  utilizing 
TS-OOK  modulation.  The  main  objective  of  our  proposed  synchronization  scheme  is  to  jointly 
determine  the  symbol  start  time  and  the  observation  window  length,  by  successively  narrowing 
down  the  integration  window  around  the  true  location  of  the  pulse.  More  specihcally,  we  consider 
that  each  packet  incorporates  a  synchronization  preamble  composed  of  a  pre-established  sequence 
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of  logical  zeros  and  ones.  During  the  preamble,  the  synchronization  block  iteratively  adapts  the 
delay  of  the  pulses  at  the  receiver  by  means  of  a  VCD  line,  which  could  be  implemented  for  example 
by  means  of  a  graphene -based  plasmonic  waveguide  [73].  The  working  principle  of  this  VCD  is 
based  on  the  property  that  the  propagation  speed  of  a  SPP  wave  on  a  graphene  waveguide  can 
be  dynamically  changed  by  electrically  modulating  the  Fermi  energy  of  the  graphene  layer  [75]. 
This  result  is  leveraged  in  our  proposed  design  to  change  the  signal  propagation  delay  between  the 
antenna  and  the  symbol  detector.  In  our  design,  we  consider  a  CTMA  detector,  which  outperforms 
energy-based  detectors  in  terms  of  sensitivity  [74]  (see  Sec.  6.2.2). 

The  proposed  synchronization  scheme,  illustrated  in  Fig.  18,  can  be  implemented  with  an  array 
of  VCD  lines  and  CTMA  detectors  and  a  single  control  logic  circuit.  Next,  we  describe  the  iterative 
procedure  by  which  the  symbol  start  time  and  observation  window  length  are  jointly  determined. 


SYNC 

signal  DATA 


Fig.  18:  Synchronization  block  diagram. 


Symbol  Start  Time  and  Observation  Window  Length  Let  m  denote  the  number  of  channels 
or  branches  in  the  synchronization  block  shown  in  Fig.  18.  Each  branch  is  composed  of  a  VCD 
line  and  a  CTMA  symbol  detector,  whose  delay  and  observation  window  length,  respectively,  are 
dynamically  adapted  through  the  following  iterative  process. 

Iteration  1:  In  the  first  iteration,  the  delay  introduced  by  each  VCD  line  is  T^/m  and  the 
CTMA  observation  window  length  is  set  to  T^/m.  Together,  all  the  CTMA  detectors  look  for  the 
pulse  over  the  entire  symbol  time  Ts.  Based  on  this,  the  output  of  the  first  CTMA  detector  is  equal  to 
one  only  when  the  pulse  is  located  between  [(m  -  IjT, /m,  T 5].  Similarly,  the  output  of  the  second 
CTMA  detector  is  equal  to  one  only  when  the  pulse  is  located  between  [(m  -  2)Tslm,  (m  -  l)r^/m], 
and  so  on.  The  output  of  the  last  branch,  m,  is  only  equal  to  one  when  the  pulse  is  located  between 
[0,  Tslm\.  We  then  define  block!  as  the  block  where  the  pulse  is  present. 

Iteration  2:  In  this  iteration,  each  CTMA  detector  utilizes  an  observation  window  of 
duration  T^/m^.  The  first  VCD  line  introduces  a  delay  equal  to  (block!  -  !)Tslm  +  Ts  I  m^.  The 
remaining  VCD  lines  introduce  a  delay  equal  to  Thus,  the  total  delay  in  the  second  line 

is  equal  to  (block!  -  !)Tslm  -t  2Ts/m^.  Similarly,  the  total  delay  in  the  third  line  is  equal  to 
(block!  -  !)Tslm  +  ^TsjrrP',  and  so  on.  We  then  define  block2  as  the  block  where  the  pulse  is 
found. 

Iteration  3:  In  this  iteration,  each  CTMA  detector  utilizes  an  observation  window  of  duration 
Ts/m^ .  The  first  VCD  line  introduces  a  delay  equal  to  (block!  -  !)Ts/m  +  (block2-!)Ts/m^  +  Ts/m^. 
The  remaining  VCD  lines  introduce  a  delay  equal  to  Ty/m^.  Thus,  the  total  delay  in  the  second  line 
is  equal  to  (block!  -  !)Tslm  +  (block2  -  !)Tslm^  -t  2Ts/m^.  Similarly,  the  total  delay  in  the  third 
branch  is  equal  to  (block!  -  !)Tslm  +  (block2  -  !)Tslm^  -l-  3Tslm^,  and  so  on.  We  then  define 
blocks  as  the  block  where  the  pulse  is  found. 
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Having  explained  the  logic  behind  the  first  three  iterations,  we  can  write  a  generalized  formula 
for  the  estimated  symbol  start  time  Td  after  the  synchronization  preamble: 


A  1  rj-i 

Td  =  'Viblocki  -  1)^  +  (blockx)^,  (56) 

■f-r*  m‘ 


where  m  is  the  total  number  of  branches  in  the  synchronization  block,  x  is  the  number  of 
synchronization  iterations  and  Ts  is  the  symbol  time.  The  maximum  number  of  branches  m  is  limited 
by  the  hardware  constraints  of  nano-devices,  and  the  number  of  iterations  x  depends  on  the  length 
of  the  preamble  bits  sent  by  the  transmitter.  Both  are  considered  parameter  values  in  our  analysis. 
Note  there  is  no  synchronization  adjustment  if  the  preamble  bits  are  equal  to  zero. 

Similarly,  the  observation  window  length  T^jin  of  the  symbol  detector  after  the  synchronization 
preamble  is  given  by 


T  —  ^ 

^  win  ~  T? 


and  can  also  be  written  as  a  multiple  of  the  pulse  duration  Tp 


(57) 


nT, 


p- 


By  combining  (57)  and  (58),  we  can  write  n  as 


n  = 


(58) 


(59) 


where  [i  is  the  spreading  factor  in  TS-OOK. 

From  Fig.  18,  once  the  synchronization  phase  is  completed,  we  close  the  DATA  out  switch  and 
disconnect  the  feedback  to  channels  two  to  m.  Hence,  we  can  use  the  same  block  to  capture  the  data 
by  using  only  the  first  branch,  with  the  VCD  delay  set  to  Td  and  the  observation  window  set  to  Tu,in- 

Impact  of  Clock  Skew  The  clock  skew,  or  difference  between  the  transmitter  and  receiver  clock 
frequencies,  introduces  further  synchronization  challenges.  In  particular,  a  mismatch  in  the  symbol 
rates  between  the  two  nodes  might  progressively  move  the  symbols  outside  of  the  detector  observation 
window.  This  sets  an  upper  bound  on  the  number  of  data  bits  that  can  be  transmitted  before  a  new 
synchronization  phase  is  needed.  In  this  section,  we  investigate  such  bound  for  different  values 
of  the  clock  skew.  In  our  analysis,  we  define  the  clock  skew  from  the  receiver’s  perspective.  In 
particular,  we  consider  that  the  symbol  time  Ts  varies  for  each  sample  n  at  a  rate  n 

Ts{n\  =  Ts{n-\]{\+n),  (60) 


where  ju  can  be  either  positive  or  negative.  When  ju  is  positive,  the  symbol  period  appears  to  increase 
with  time,  whereas  when  jj.  is  negative,  it  appears  to  decrease. 

Without  loss  of  generality,  we  consider  that,  after  the  initial  synchronization  phase,  the  position 
p  of  the  pulse  within  the  observation  window  is  characterized  with  a  uniform  random  variable  with 
probability  density  function  (PDF)  fp  given  by 


fpip)  = 


0, 


if  0  <  ;?  <  ^, 
otherwise. 


(61) 
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T  T 

Therefore,  the  expected  pulse  position  is  given  by  or  equivalently,  The  maximum  number 
of  bits  Hbits  to  be  transmitted  before  requiring  a  new  synchronization  phase  is  then  computed  as 
follows.  When  jj.  is  positive,  re- synchronization  should  happen  before  this  condition  is  met 


Ts[n\  >  r^[0]  +Tuji„l2 


(62) 


where  /r  is  the  clock  variation  factor,  is  the  symbol  time  and  Tp  is  the  pulse  time.  By  combining  (62), 
(60)  and  (59),  nbus  can  be  obtained  as 


r,[0]-(i  +  /r)'’^-  =(r,[0]  + W2), 
_  logio(l  + 
logio(l+ju) 


(63) 

(64) 


We  can  see  that,  as  the  value  of  jx  increases,  the  number  of  bits  after  which  re-synchronization 
should  take  place  decreases,  as  expected.  The  same  procedure  can  be  followed  to  obtain  the  value  of 
Hbits  for  /r  <  0,  reaching  the  same  equation  for  Hbus  in  (64),  but  with  |/i|  instead  of  yu. 


6.2.2  Continuous-Time  Moving-Average  Detector 

In  this  section,  we  describe  and  analytically  model  the  performance  of  the  CTMA  symbol  detector 
scheme,  with  a  special  emphasis  on  the  impact  of  observation  window  on  the  successful  symbol 
detection  probability. 

Detector  Overview  Traditional  energy-based  symbol  detector  schemes  for  pulse-based  communi¬ 
cation  integrate  the  received  signal  over  the  entire  symbol  time  or,  following  our  definition,  over  the 
observation  window  Twin-  The  output  is  then  compared  to  a  pre-defined  threshold  and  a  decision 
on  the  received  symbol  is  made.  However,  as  mentioned  before,  the  signals  radiated  by  individual 
nano-devices  have  a  very  low  power  and  are  affected  by  the  THz-band  channel  path-loss.  This 
results  into  very  high  Symbol  Error  Rate  (SERs)  in  realistic  scenarios. 

To  overcome  such  limitation,  a  CTMA  symbol  detector  for  THz  pulses  was  first  proposed 
in  [74].  A  revised  block  diagram  for  the  CTMA  symbol  detector  is  shown  in  Eig.  19. 

SIGNAL  SYMBOL 


Start  Stop 


Fig.  19:  CTMA  symbol  detector  block  diagram. 

The  fundamental  idea  behind  this  symbol  detection  scheme  can  be  summarized  as  follows. 
When  the  detector  integrates  the  pulse  over  a  window  Twin  greater  than  the  pulse  time  Tp,  the  noise 
contribution  increases  with  Twin-  This  effect  drops  the  performance  of  the  receiver  in  terms  of  SER 
(Eig.  20(top)).  To  mitigate  the  impact  of  noise,  we  could  reduce  the  integration  time  down  to  the 
pulse  duration,  Tp.  However,  since  the  observation  window  T^m  is  generally  much  larger  than  Tp, 
and  we  cannot  know  the  exact  position  of  the  pulse  within  T^m,  we  would  have  to  implement  a  total 
of  TwinjTp  integrators  in  parallel  and  shifted  in  time  to  guarantee  that  the  symbol  is  not  missed.  This 
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implementation  is  certainly  impractical  for  low  complexity  devices.  Additionally,  the  pulse  may  not 
perfectly  fit  in  the  integration  time  window  (Fig.  20  (center)).  Alternatively,  we  could  just  use  a 
single  integrator  and  move  it  in  time  at  tiny  steps  At  ^  0.  This  is  the  concept  of  a  CTMA  symbol 
detector  (Fig.  20(bottom)). 


Fig.  20:  Window  length  bigger  than  Tp  (top);  Window  length  equal  to  Tp,  but  pulse  does  not  fit  in  a  single 
integration  window  (center);  Window  length  equal  to  Tp  and  pulse  fits  in  a  single  integration  window  (bottom). 


We  define  the  CTMA  detector  as  the  linear  time-invariant  system  with  input  to  output  relation 
given  by 

t 

x{t)  =  J  v(t)dt,  (65) 

t—Tp 

where  x  stands  for  its  output,  t  refers  to  time,  Tp  is  the  pulse  time  and  v  is  the  input  signal  to  CTMA. 
As  discussed  in  [74],  the  CTMA  behaves  as  a  low-pass  filter  and  could  be  implemented  with  a 
discrete  analog  components.  At  each  instant  of  time,  the  output  of  the  CTMA,  x  is  compared  to  a 
fixed  threshold  Vth,  and  the  maximum  value  is  held  until  the  output  is  read  and  reset.  The  detector 
then  decides  whether  the  n-th  symbol  is  one  or  zero  according  to: 


fl,  if  y  =  m2iX  {x{t))  >  Vth, 
1 0,  otherwise. 


(66) 


for  t  6  [0,  Tu,in\-  Next,  we  analytically  investigate  the  successful  symbol  detection  probability  for 
both  zeros  and  ones. 

Detection  of  Zeros  In  TS-OOK,  zero  is  transmitted  as  silence  and,  thus,  only  noise  is  detected 
at  the  receiver  side.  For  this,  we  consider  the  detected  signal  to  be  composed  of  N  identically 
distributed  random  variables.  We  define  X  =  {Xi,  X2,  ...,X7v}  as  the  vector  that  contains  the  N 
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random  variables  and  we  use  Xi  to  refer  to  the  random  variable  that  represents  the  value  of  the  input 
signal  X  at  time  i/.  N  depends  on  the  observation  window  length,  Tu;in,  the  integration  window 
length  of  the  CTMA,  which  we  consider  to  be  equal  to  Tp,  and  the  time  step  of  the  CTMA  detector, 
z.  From  (59),  N  can  be  written  as 

A  =  ^^  +  l.  (67) 

z 

Strictly  speaking,  z  ^  0,  as  the  CTMA  is  an  analog  detector  and,  thus,  N  ^  oo.  Then,  from  (66), 
the  output  of  the  symbol  detector  is  given  by  T  =  max  {X}.  Therefore,  the  PDF  fy^no  of  the  output 
can  be  obtained  as  a  function  of  the  PDFs  for  the  individual  variables  Xi.  In  particular,  it  follows  that 

fY.no  (y)  =  -^(Fno  (y)f  =  NF„o  fno  iv)  ,  (68) 

dy 

where  fno  and  Fno  stand  for  the  PDF  and  the  Cumulative  Distribution  Function  (CDF)  for  each  one 
of  the  individual  random  variables.  From  [76],  we  consider  that  each  individual  variable  follows  a 
chi-square  distribution  and,  thus. 


iv-2)l2^-iyl2) 

/.»(■<)= 

where  y  =  2X /Nq  is  the  normalized  signal  with  respect  to  the  noise  two-sided  power  spectral  density 
Nq.  V  =  2TpB  is  the  number  of  degrees  of  freedom,  where  Tp  and  B  stand  for  the  pulse  duration  and 
its  equivalent  bandwidth,  and  F  refers  to  the  Gamma  function.  Because  we  are  using  realistic  THz 
pulses,  which  are  further  affected  by  the  THz -band  channel  and  the  CTMA  hlter,  the  product  TpB  is 
generally  larger  than  one  and,  thus,  there  is  no  closed  form  solution  for  (68). 

Detection  of  Ones  In  TS-OOK,  a  logical  one  is  transmitted  as  a  one-hundred-femtosecond-long 
pulse.  Following  a  similar  methodology  as  before,  we  consider  that  the  received  signal  is  contributed 
by  two  kinds  of  random  variables,  some  representing  the  pulse  and  some  presenting  noise.  Hence, 
the  output  can  be  now  represented  as  T  =  max  {X}  =  max  {X^,  X5},  where  X and  X5  are  vectors 
that  contain  A„  noise  random  variables  and  Ns  signal  random  variables.  As  discussed  in  [74], 
provided  that  the  signal  power  is  kept  higher  than  the  noise  power,  Y  «  max  {X5}.  Therefore,  the 
PDF  of  the  output  fy^si  can  be  now  written  as 


fY.si{y)  =  NsFsi{yf^-^  fsiiy),  (70) 

where  now  each  variable  X,-  is  modeled  as  non-central  chi-square  distribution  with 

h(u)  =  \[lf  (71) 

where  A  =  =  v  SNR,  v  =  2TpB,  and  /„  stands  for  the  n-th  order  modified  Bessel  Function  of  the 

hrst  kind.  As  before,  there  is  no  closed-form  expression  for  (70). 

Symbol  Error  Rate  To  analyze  the  performance  of  the  symbol  detector  as  a  function  of  the 
observation  window  length  T^in,  we  are  interested  in  computing  the  SER.  This  is  directly  obtained  as 

SER  =  Pe\s=0Ps=0  +  Pe\s=lPs=h  (72) 
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where  po  and  pi  refer  to  the  probability  to  transmit  0  and  1,  respectively,  and 

CO 

Pe\s=o  =  /  fY,no  (?/)  dy,  (73) 

^Vth 

Nq 

2V,h 

No 

Pe\s=i  =  j  fY,si  (y)  dy,  (74) 

0 

where  fy^no  and  fy^si  are  given  by  (68)  and  (70),  respectively,  and  VthINo  stands  for  the  normalized 
signal  power  for  the  detection  of  logic  0  and  logic  1  and  No  stands  for  the  noise  level.  The  threshold 
value  is  calculated  by  finding  out  the  intersection  of  logic  0  and  logic  1. 

6.2.3  Throughput  Analysis 

From  the  discussion  in  Sec.  6.2.2,  we  are  interested  in  reducing  the  observation  window  length  T^jin 
to  reduce  the  value  of  N  in  (68)  and  the  SER  given  by  (72).  In  Sec.  6.2.1,  we  showed  that  with  the 
proposed  synchronization  block.  Twin  can  be  reduced  by  increasing  the  number  of  synchronization 
bits  X  in  (57).  Transmitting  longer  preambles  can  result  in  a  reduction  of  the  throughput,  but  at  the 
same,  given  that  the  SER  is  reduced,  it  is  more  likely  that  the  transmitted  bits  are  properly  received, 
which  improves  the  throughput.  Therefore,  there  is  a  compromise  between  different  parameter 
values  on  the  synchronization  and  symbol  detection  block. 

In  this  section,  we  aim  at  identifying  the  best  values  for  such  parameters,  by  defining  fhe 
maximization  of  the  throughput  as  an  optimization  problem.  In  particular,  the  problem  is  formulated 
as: 


Given: 

Eind: 

Maximize: 


Subject  to: 


where 


SNR,fi 


in,  X,  nyiis 


T delay 


loglo(l  +  2i 


in  <  (5,  nhits  —  ,  I  ’ 

logio(l  +  |//|) 

SER  =  ^(m,  ;c,  SNR)', 


FER  =  A  {SER,  nuts,  x);  y  =  l/FER; 

P delay  ~  ~  1)  i^bits  P^ , 


(75) 


•  The  first  constraint  guarantees  that  the  final  observation  window  length  T^in  is  shorter  than 
the  symbol  duration  Ts,  as  given  by  (59). 

•  The  second  constraint  guarantees  that  there  is  no  need  for  re-synchronization  within  the 
transmission  of  one  frame.  This  depends  on  the  clock  skew  p  in  (64). 
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•  'F  in  the  third  constraint  represents  the  SER  as  a  function  of  m,  x  and  SNR,  and  is  given 
by  (72). 

•  A  in  the  fourth  constraint  represents  the  Frame  Error  Rate  (FER)  and  is  a  function  of  the 
SER  and  the  total  number  of  bits  per  frame,  including  the  synchronization  preamble.  If  we 
consider  that  a  simple  error  detection  scheme  is  used,  without  error  correction  capabilities, 

A  =  (1  -  (76) 

where  P  is  the  total  number  of  bits  in  the  preamble,  P  >  x,  as  not  all  the  bits  in  the  preamble 
might  be  ones. 

•  rj  refers  to  the  average  number  of  retransmissions  per  frame  and  T^eiay  refers  to  the  average 
delay  per  frame. 

The  optimal  values  for  m,  x  and  n„its  that  maximize  the  single  link  throughput  can  be  then 
numerically  obtained. 

6.3  Results  and  Discussions 

In  this  section,  we  hrst  test  the  functioning  of  the  proposed  synchronization  scheme  with  experimental 
THz  pulses.  Then,  we  numerically  investigate  the  impact  of  different  parameters  on  the  performance 
of  the  synchronization  and  detection  blocks  and  compute  the  achievable  throughput. 

6.3.1  Synchronization  Test  with  Experimental  Data 

To  test  the  functioning  of  the  proposed  synchronization  algorithm,  we  utilize  measured  THz  pulses 
generated  and  captured  by  means  of  a  time-domain  THz  spectroscopy  platform  available  to  the 
authors.  Such  platform  is  able  to  generate  one-hundred-femtosecond-long  pulses,  which  are  emitted 
in  free  space  and  measured  at  a  distance  by  a  detector.  The  captured  signal  is  passed  through  the 
proposed  synchronization  block  in  Fig.  18  which  iteratively  estimates  the  symbol  start  time  and, 
thus,  the  necessary  delay  to  be  introduced  by  the  first  VCD  during  the  reception  of  the  actual  data 
symbols. 

The  different  iterations  of  the  algorithm  for  m  =  .r  =  2  are  illustrated  in  Fig.  21.  In  the  first 
iteration,  CTMAi  looks  from  Ts  jl  to  Ts  and  CTMA2  looks  from  0  to  Ts  jl.  CTMA2  detects  the 
pulse  and,  thus,  Blockl=2.  In  the  second  iteration,  CTMAi  looks  from  r^/4  to  r^/2  and  CTMA2 
looks  from  0  to  TsjA.  Now  CTMAi  detects  the  pulse  and,  thus,  Block2=l.  During  the  data  phase, 
VCDi  is  tuned  to  introduce  a  delay  equal  to  0.5Ts  and  CTMAi  utilizes  an  observation  window 
length  equal  to  Tj/d.  In  Table  2,  we  summarize  the  output  for  different  number  of  pulses  in  the 
synchronization  preamble  x. 

After  the  synchronization  preamble,  the  measured  data  symbols  were  detected  by  utilizing  the 
estimated  delay  for  VCDi  and  the  observation  window  length  for  CTMAi.  The  low  noise  in  the 
measurements  resulted  in  almost  no  symbol  errors,  in  line  with  the  numerical  results  discussed  next. 
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Fig.  21:  Illustration  of  the  synchronization  algorithm  for  m  =  x  =  2:  First  iteration  (top);  Second  iteration 
(middle);  Third  iteration  (bottom). 

Table  2:  Iterations  in  the  synchronization  algorithm  (m=2) 


Preamble 

it=l 

it=2 

it=3 

it=4 

Td 

Twin 

X  =  2 

2 

2 

NA 

NA 

Ts 

Tsl4r 

X  =  3 

2 

2 

2 

NA 

Ts 

Ts/S 

x  =  4 

2 

2 

2 

1 

iSTsjie 

Ts/ie 

6.3.2  Numerical  Results 

We  consider  the  following  parameters  in  our  analysis.  The  transmitted  pulses  are  modeled  as  the 
first  time  derivative  of  a  Gaussian  pulse  with  pulse  length  Tp=l00  fs  and  pulse  energy  Ep=l  aJ. 
The  time  between  symbols  Ts  is  set  to  10  ps  =  100).  Symbols  are  considered  equiprobable,  i.e., 
Ps=o  =  Ps=i  =  0.5  in  (72). 

Impact  of  Clock  Skew  on  the  Frame  Length  In  Fig.  22,  the  number  of  bits  nhns  after  which 
re  synchronization  is  needed  given  by  (64)  is  shown  as  a  function  of  number  of  branches  in  the 
synchronization  block  m,  for  different  values  of  the  clock  skew  p.  As  expected,  for  higher  values  of 
p,  the  number  of  bits  after  which  re- synchronization  should  take  place  decreases  hence  signifying, 
re-synchronization  becomes  more  often  and  more  preamble  bits  are  needed  to  transmit  the  same 
amount  of  data  bits.  Moreover,  as  the  number  of  channel  increases,  the  number  of  bits  after  which 
re-synchronization  should  take  place  decreases.  This  is  because,  T^jin  decreases  as  we  increase  m 
and  hence,  the  CTMA  looks  to  re-synchronize  again  when  the  clock  shifts  by  a  value  of  T^in  which 
is  small  for  large  values  of  m. 

Symbol  Error  Rate  In  Fig.  23,  the  SER  given  by  (72)  is  shown  as  a  function  of  the  SNR  for  four 
different  values  of  m^.  As  the  term  increases,  the  observation  window  length  Twin  decreases  and. 
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Fig.  22:  Maximum  number  of  bits  before  resyncbronization  as  a  function  of  ,  for  different  values  of  the 
clock  skew  /i. 

thus,  the  SER  also  decreases.  This  can  also  be  understood  from  the  PDFs  for  the  detected  signal  in 
the  transmission  of  Os  and  Is,  given  by  (68)  and  (70),  respectively. 


Fig.  23:  SER  as  a  function  of  the  SNR  for  different  values  of  . 

Throughput  In  Fig.  24,  the  throughput  given  by  (75)  is  plot  as  a  function  of  the  for  different 
values  of  the  SNR.  It  can  be  seen  that  as  the  value  of  increases,  the  throughput  also  increases. 
This  result  shows  that  the  benefits  of  reducing  SER  by  shortening  the  observation  window  length 
Tuiin  largely  overcome  the  cost  introduced  by  increasing  the  size  of  the  synchronization  block,  i.e., 
increasing  m,  or  transmitting  longer  synchronization  preambles,  x. 

Given  the  constraints  of  nano-devices,  it  seems  more  adequate  to  keep  m  to  the  minimum 
number  possible,  i.e.,  m  =  2,  and  investigate  the  impact  of  x  alone.  In  Fig.  25,  the  throughput  is 
shown  as  a  function  of  the  v  for  m  =  2.  For  low  SNR  values,  the  impact  of  x  on  the  SER  is  larger 
than  the  impact  of  x  in  the  total  delay.  Therefore,  the  throughput  S  increases  rapidly  with  x  at  low 
values  of  SNR.  For  example,  for  SNR=3  dB,  we  can  triple  the  throughput  by  adding  only  6  bits 
to  the  synchronization  preamble.  For  higher  values  of  SNR,  the  impact  of  x  on  the  SER  is  less 
significant  and  hence  the  throughput  does  not  drastically  change  with  x. 

Finally,  the  impact  of  the  number  of  bits  per  frame  on  the  throughput  S  is  illustrated  by  means 
of  Fig.  26.  For  low  values  of  SNR,  lesser  of  bits  should  be  sent  in  a  single  packet,  whereas  for  large 
values  of  SNR,  the  trend  reverses  and  the  throughput  remains  the  same  for  the  different  numbers  of 
bits. 

Optimal  Synchronization  and  Detection  Parameters  In  light  of  the  illustrated  results,  combined 
with  exhaustive  search,  the  values  of  m,  n,  SNR  and  bits  that  are  required  to  maximize  the  throughput 
can  be  summarized  as  follows.  For  SNR  above  7  dB,  the  throughput  remains  almost  constant 
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Fig.  24:  Throughout  as  a  function  of  the  SNR  for  different  values  of  (frame  length=1000  bits). 


Fig.  25:  Throughout  as  a  function  of  the  SNR  for  different  values  of  v  (m=2,  frame  length=1000  bits). 

for  different  values  of  x,  m  and  bits.  In  that  ease,  redueing  m  and  x  allows  the  transmission 
of  longer  frames,  if  needed,  and  simplifies  the  synchronization  block  hardware.  Otherwise,  the 
maximum  throughput  is  numerically  obtained  for  x  =  6,  m  =  3  and  bits  =  5000.  For  low  SNR 
values,  significant  improvement  is  achieved  by  increasing  the  number  of  bits  in  the  synchronization 
preamble  x.  Maximum  throughput  is  achieved  again  for  x  =  6  and  m  =  3,  now  with  shorter  frames. 

6.3.3  Discussion 

We  have  a  proposed  a  joint  physical-layer  synchronization  and  symbol  detection  scheme  for 
pulse-based  THz-band  communication.  The  main  objective  of  the  proposed  scheme  is  to  jointly 
determine  the  symbol  start  time  and  the  observation  window  length,  by  successively  narrowing 
down  the  integration  window  around  the  true  location  of  the  pulse.  The  proposed  solution  is 


Fig.  26:  Throughout  as  a  function  of  the  SNR  for  different  values  of  bits  {m=3,  x=6). 
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fully  analog  and  can  be  implemented  with  a  combination  of  voltage-controlled  delay  lines  and 
continuous-time  moving-average  symbol  detectors  with  tunable  observation  window.  We  have 
analytically  investigated  the  synchronization  preamble  length  and  the  maximum  number  of  bits  to 
be  transmitted  before  requiring  re-synchronization.  We  have  also  obtained  analytical  expressions 
for  symbol  error  rate  as  a  function  of  the  observation  window  length.  We  have  investigated  the 
compromise  between  synchronization  preamble  length  and  symbol  detection  observation  window 
length,  and  investigated  their  joint  impact  on  the  achievable  throughput.  Finally,  we  have  tested  the 
algorithm  with  experimental  THz  pulses  and  provided  extensive  simulation  results  to  analyze  the 
impact  of  different  parameters.  The  results  show  how,  with  the  proposed  synchronization  algorithm, 
tight  synchronization  and  low  symbol  error  rates  are  possible  just  with  very  short  preambles,  less 
than  bits  long. 
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7  PACKET  SIZE  OPTIMIZATION 
7.1  Motivation 

One  of  the  key  applications  of  THz  communications  is  in  the  field  of  nanonetworks  and  WNSNs  [29, 
34].  Indeed,  the  very  small  size  of  graphene-based  plasmonic  nano-transceivers  and  nano-antennas 
(Sec.  4)  enables  their  integration  with  miniature  sensors  for  physical,  chemical  and  biological  sensing. 
On  the  one  hand,  the  THz  band  provides  nanosensors  with  a  very  large  bandwidth  [23].  However, 
on  the  other  hand,  the  very  high  propagation  loss  at  such  frequencies  combined  with  the  very  limited 
transmission  power  of  nanosensors  result  in  very  short  communication  distances.  Moreover,  the 
limited  capacity  of  nano-batteries,  which  requires  the  use  of  time-consuming  energy-harvesting 
procedures  [35,36],  and  the  limited  computational  capabilities  of  nanosensors,  affect  the  throughput 
of  WNSNs.  All  these  interdependencies  motivate  the  joint  analysis  of  the  nano-device  capabilities, 
the  THz  band  peculiarities  and  their  impact  on  the  achievable  throughput. 

There  have  been  many  cross-layer  studies  on  packet  size  optimization  in  wireless  commu¬ 
nication  networks  for  a  myriad  of  environments,  including  terrestrial,  underwater,  underground, 
and  intra-body  sensor  networks  [77-80].  However,  all  these  works  cannot  be  directly  applied  for 
energy-harvesting  networks,  in  which  the  energy  fluctuates  with  time  instead  of  monotonically 
decreasing.  In  this  direction,  several  energy  consumption  optimization  problems  for  wireless 
networks  with  energy  harvesting  nodes  have  been  proposed  over  the  recent  years  [81-83].  All 
these  works  are  mainly  focused  on  optimizing  the  utilization  of  the  harvested  energy  following 
a  general  approach  to  And  the  trade-off  between  the  consumed  energy  and  the  achieved  quality 
of  service.  While  the  aforementioned  studies  are  applicable  to  general  wireless  communication 
networks,  in  [84],  a  study  is  performed  for  the  specific  case  of  communication  in  nanonetworks. 
However,  the  impact  of  such  energy  management  policies  on  the  achievable  throughput  at  the  link 
layer  is  not  analyzed. 

As  part  of  this  effort,  we  have  addressed  the  throughput  maximization  problem  in  WNSNs, 
by  taking  into  account  the  existing  device  and  communication  interdependencies.  In  particular,  a 
link  throughput  optimization  problem  is  defined,  and  the  optimal  data  packet  size  which  maximizes 
the  link  efficiency  is  derived  by  capturing  the  power,  energy  and  computational  constraints  of 
nanosensors;  the  very  high  path-loss  and  very  large  bandwidth  of  the  THz-band  channel;  the 
possibility  to  communicate  by  transmitting  one-hundred-femtosecond-long  pulses,  which  can 
virtually  create  parallel  orthogonal  channels  between  nanosensors  [63];  and  three  different  error 
control  strategies  tailored  to  WNSNs,  namely.  Automatic  Repeat  reQuest  (ARQ),  Forward  Error 
Correction  (EEC)  and  novel  Error  Prevention  Codes  (EPCs),  which  have  been  designed  with  the 
WNSNs  peculiarities  in  mind  [85].  Both  the  energy  harvesting  limits  and  the  successful  packet 
transmission  time  are  defined  as  the  optimization  problem  constraints,  and  the  optimal  solution  is 
derived  by  using  a  bisection  method.  Our  results  show  that  the  link  efficiency  quickly  decreases  when 
considering  the  energy  constraints  compared  to  the  scenario  that  there  is  no  energy  shortage.  The 
decrease  depends  on  various  parameters  including  the  error-control  technique,  the  communication 
distance,  and  the  harvesting  capability  of  the  nanosensors.  Similarly,  the  packet  size  quickly 
decreases  with  the  transmission  distance,  approaching  several  hundreds  bits  for  distances  beyond  a 
few  millimeters. 

Next,  in  Sec.  7.2.1,  we  define  the  system  model  for  WNSNs  and  discuss  the  throughput 
optimization  problem  while  defining  different  constraints  for  various  error-control  methods.  Sec.  7.2.2 
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contains  the  approach  to  solve  the  optimization  problem,  and  covers  the  related  algorithm  to  find  the 
optimal  solution  for  our  problem.  Finally,  in  Sec.  7.3.1,  we  numerically  study  the  optimal  packet 
size  for  different  error-control  techniques  and  compare  the  performance  of  them  under  different 
conditions. 


7.2  Methods,  Assumptions  and  Procedures 

7.2.1  System  Model  and  Optimization  Problem  Formulation 


We  consider  a  single  communication  link  between  two  nanosensors.  Nanosensors  communicate 
with  each  other  using  TS-OOK  [63].  As  previously  discussed,  such  mechanism  can  effectively 
provide  nanomachines  with  orthogonal  communication  challenges,  thus  minimizing  the  potential 
multi-user  interference.  The  transmitter  node  is  a  nanomachine  with  the  capability  of  harvesting 
energy  by  means  of  piezoelectrie  nano-generators  [35],  whieh  converts  kinetic  energy  into  eleetrieity 
by  exploiting  nanowires.  At  the  receiver  node  it  is  assumed  that  always  enough  energy  is  available 
to  receive  the  packets  successfully,  and  the  receiver  has  enough  amount  of  memory  to  buffer  fhe 
reeeived  data.  Based  on  these  assumption,  first  a  throughput  optimization  problem  with  the  objective 
of  maximizing  the  link  throughput  between  a  pair  of  transceivers  in  a  WNSN  is  defined  in  this 
section,  considering  different  constraints  for  different  error-control  techniques.  Then  different 
eonstraints  are  tailored  in  details  for  different  error-eontrol  teehniques  respeetively.  We  start  with  the 
definition  of  throughput  in  a  WNSN  link  as  the  rate  of  sueeessful  message  transmission  measured  in 
bits  per  second  (bps),  which  can  be  given  by: 

^  Useful  Data  Length 

Successful  Transmission  Time 

Without  loss  of  generality,  instead  we  ean  optimize  the  link  utilization  efficiency  rj  for  a  given 
transmission  rate  r  which  can  be  further  defined  as: 


Ldata 
Nret  •  Ttx 


(78) 


where  Ldata  is  the  useful  data  length  in  bits,  Nret  represents  the  expeeted  number  of  retransmissions 
needed  for  the  packet  to  be  received  and  processed  successfully  at  the  receiver  node  according  to 
the  chosen  error-control  method.  Ttx  is  the  total  time  required  to  accomplish  a  complete  packet 
transmission  including  the  time  required  to  harvest  enough  energy  for  transmission,  and  will  be 
defined  later  in  See.  7.2.1  as  one  of  the  eonstraints  of  the  optimization  problem.  Henee,  for  a  chosen 
error-control  technique,  to  find  optimal  data  length  which  maximizes  the  link  efficiency,  the 
optimization  problem  can  be  formulated  as  follows. 

Channel  Efficiency  Maximization  Problem  for  WNSNs  We  define  next  the  optimization  prob¬ 
lem  with  the  objeetive  funetion  of  maximizing  the  ehannel  efficieney  in  WNSNs,  and  the  eonstraints 
whieh  will  be  defined  later  in  See.  7.2. 1  through  7.2. 1  in  details  aceording  to  the  error-eontrol  method. 


Optimization  Problem  [PI] 
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Given  : 

r,  (f,  Nq,  P%  f,  d,  6^,  T^,  Ef^,  Aharv 

Find  : 

L* 

data 

Maximize  : 

1  Ldata 

^  r  Nret  •  Ttx 

Subject  to  : 

Ldata  ^  0 

(79) 

1 

Nret  —  ^ 

1  -Pe 

(80) 

Pe  =  (BER,  Ldata,  L^) 

(81) 

BER  =  cp^(^^,p^) 

^PL-Nq  ’ 

(82) 

Ttx  = 

(83) 

rpQ  _  n&  ( j  j  G  yfi  j,G  j,  'V 

^tx~^  \‘-‘data,‘-‘  ,  X  proc^  ^  tjo' ^  P''°P' 

(84) 

'^tx-harv  ~  (Ldata,  L  ,  p  Aharv) 

(85) 

Efx  ^  '^harv  '  Lf^-harv 

(86) 

where: 

•  Pe  =  iJ/^(BER,  Ldata,  L^)  IS  the  Packet  Error  Rate  (PER),  which  is  a  function  of  the  Bit  Error 
Rate  (BER),  data  length  Ldata,  and  the  length  of  the  redundant  bits  L^,  which  depends  on  the 
error-control  scheme  S. 

•  BER  =  is  a  function  of  energy  required  to  transmit  a  bit  E^^^,  the  path  loss 

between  the  transmitter  and  receiver  nodes  PL,  and  the  noise  spectral  density  Nq.  It  also 
depends  on  the  pulse  probability  pf  that  itself  varies  according  to  the  chosen  error-control 
method  S. 

•  PL(  f,  d)  is  the  path  loss  between  the  transmitter  and  receiver  nodes  and  depends  on  the 
transmission  frequency  /  and  the  distance  between  the  nodes  d. 

•  ~  L^,Tproc^T'fla^'^prop)  is  the  packet  round-trip  time  or  the  total  time  needed  for 

the  packet  to  be  transmitted  and  the  acknowledgment  to  be  received.  and  are  the 
processing  time  for  a  complete  packet  transmission  and  the  time-out  before  retransmission 
respectively,  and  will  be  later  defined  in  Sec.  7.2.1  through  7.2.1  based  on  chosen  &.  Tprop 
is  the  propagation  time  and  is  relative  to  the  distance  between  the  transmitter  and  receiver 
nanomachines  and  it  can  be  rewritten  as  djc  where  c  is  the  speed  of  light. 

•  (Ldata,  L^,  Ef^,p^,  Aharv)  Is  the  time  to  harvest  enough  energy  for  a  complete 
transmission  including  processing  and  transmitting  the  data. 
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•  Ef^  is  the  energy  eonsumption  needed  to  transmit  a  paeket  of  length  L^ata-  Note  that  this 
energy  eontains  both  the  required  energy  to  proeess  data  whieh  depends  on  seleeted  fi,  and 
the  required  energy  for  transmitting  the  data  whieh  depends  on  Ldata,  E^ ,  and  pf . 

•  ^harv  is  the  rate  at  whieh  the  nanosensor  transmitter  mote  is  able  to  harvest  energy  in  J/s. 

Here  we  assume  that  the  nano-transmitter  always  transmits  with  the  maximum  available 
energy,  henee  the  inequality  in  (86)  reduees  to  the  following  equality: 

Efx  =  ^harv  •  T^x-harv’’ 

Therefore  the  funetion  r^,  ean  be  dehned  as  follows: 

E^ 

rfi  ^ 

^harv 

Henee  the  two  eonstraints  (85)  and  (86)  ean  be  merged  into  one  eonstraint  as  in  (88).  Moreover,  for 
the  BER,  the  funetion  does  not  depend  on  the  optimization  variable  Ldata,  and  depends  on  the 
physieal  layer  parameters.  In  this  paper  we  use  the  derived  values  for  BER  in  [85],  as  given  values 
of  the  proposed  optimization  problem. 

As  it  ean  be  seen,  the  optimization  problem  dehned  in  [PI]  is  a  general  problem  with  funetions 
,  9^,  and  Ep^  whieh  have  to  be  dehned  for  three  different  error-eontrol  teehniques,  namely  ARQ, 
FEC,  and  EPC,  whieh  are  addressed  as  follows. 

Automatic  Repeat  reQuest  (ARQ)  Constraints  The  packet  error  rate  in  ARQ  is  dehned  as 
follows: 

=  1  -  (1  -  BER^^^y,  (89) 

where  I  =  Ldata  +  Lcrc,  and  Lcrc  is  the  length  of  Cyclic  Redundancy  Check  (CRC)  used  for  error 
deteetion.  Moreover,  the  packet  round-trip  time  is  given  by  the  following  equation: 


rpARQ 
^  tx 


rpARQ  rp 

tx,data  CRC 
ARQ  ARQ, 


pARQ 


+  Ps,dataPs,ack(^'^P>-^P  +  '^CRC  + 
+ (1 _ 

'  s,data‘^ s,ackt  tjo  ’ 


(90) 


where  and  are  the  data  and  aeknowledgment  transmission  times  and  are  given  by  I  jr 

tx,data  tx,ack  °  o  j  / 

and  Lackjc  respeetively,  where  Lack  is  the  aeknowledgment  paeket  length.  Tcrc  refers  to  the  delay 
eaused  by  eomputational  proeess  of  the  CRC,  and  is  given  by  {Ldata  •  Tdk)  where  Tdk  is  the  inverse  of 
the  nanomaehine’s  eloek  frequeney.  p'^data  P^afk  acknowledgment  paeket  sueeess 

rate  respeetively  and  ean  be  obtained  from  (89)  with  I  =  Ldata  and  I  =  Lack  respeetively.  Finally 
is  dehned  as  follows: 

tjo 


=  i.i(2v„, + w + )•  (91) 

whieh  is  the  propagation  delay  to  transmit  the  data  and  reeeive  the  aeknowledgment,  plus  the  time 
the  reeeiver  takes  to  proeess  the  CRC  and  transmit  the  aeknowledgement  paeket,  plus  a  ten  pereent 
margin  time.  The  energy  required  for  transmission  proeess  in  ARQ  ean  be  dehned  as  follows  [85]: 
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(92) 


^CRC  + 

tx  tx,aata  '-jv'-  s,aata^  s,ack  rx,ack^ 

where  and  ,  refer  to  the  energy  required  to  transmit  the  data  paeket  and  reeeive  the 

aeknowledgement  paeket,  and  are  given  by  {I  •  p^  •  Ef'J)  and  •  P^^fk  '  •  Ej^^J)  respeetively. 

Ecrc  stands  for  the  eonsumed  energy  eaused  by  eomputational  proeess  of  the  CRC,  and  is  given  by 
(Lcrc  ■  LdataiEshift  +  Ehoid))  where  Eshi  ft  and  Ehoid  are  the  energy  eonsumed  to  shift  and  hold  a 
registry  value  in  a  shift  register.  Note  that  the  value  of  Ecrc  is  defined  based  on  the  assumption 
that  a  CRC  is  implemented  by  exploiting  shift  registers  and  XOR  logie  gates  as  deseribed  in  [85]. 
Therefore,  the  equivalent  optimization  problem  for  the  ARQ  mode,  is  given  by  [PI],  while  replaeing 
the  eonstraints  (81)  and  (84)  by  (89)  and  (90),  respeetively,  and  using  eonstraint  (88)  instead  of  (85) 
and  (86),  while  the  definition  of  Ef^  is  given  by  (92). 

Forward  Error  Correction  (FEC)  Constraints  For  the  FEC,  depends  on  the  Block  Error 
Rate  (BEER),  and  is  defined  as  follows: 


r 


where  n  is  the  number  of  bloeks  per  data  paeket  payload,  and  BLER^^^  is  given  as  follows: 


BLER^^^  =  ^  [M(BER)^(1  -  BERf-j, 


j-t+\ 


(93) 


(94) 


where  k  refers  to  the  bloek  size  and  t  is  the  error  eorreetion  eapability  of  the  eode.  Moreover,  the 
packet  round-trip  time  for  EEC  ean  be  defined  as  follows: 


't'FEC  _  't'FEC  I  't'FEC  I  EEC  /'t<  .  't<FEC  \ 

^tx  -  ^tx4ata  ^code  Ps^atoS^P^op  ^decode) 


+  (^-plStaFtlo 


(95) 


where  T^^ata  is  the  data  transmission  time  and  is  given  by  I  jr  where  I  is  the  total  length  of  the 
transmitted  data  and  is  equal  to  {Ldata  +  and  is  the  length  of  the  redundant  bits  added 

for  error  eorreetion.  and  refer  to  lateney  eaused  by  eoding  and  deeoding  proeesses  of 

the  data  respeetively  and  are  given  as  follows  [85]: 


jFEC 
^  code 


=  2nTcik 


(96) 


jFEC 
^ decode 


=  {k  +  l)nTcik, 


(97) 


P^sdata  paeket  sueeess  rate  and  ean  be  obtained  from  (93).  Einally,  is  given  as 

follows: 


rFEC 

tjo 


1.1  (2r, 


prop 


,  jEEC  'I 

decode^ 


The  funetion  Efl  for  EEC  mode  is  defined  as  follows: 


^tx 


pFEC  _  pFEC  ,  pEEC 
^tx  ^tx,data  ^code  ’ 


(98) 


(99) 


Approved  for  Public  Release;  Distribution  Unlimited 

47 


where  refers  to  the  energy  required  to  transmit  the  data  packet,  and  is  given  by  {I  ■ 

and  stands  for  the  consumed  energy  caused  by  computational  process  of  coding  the  data  in 
transmitter,  and  is  given  by  n  ■  k{Eioad  +  Ehoid)  where  Eioad  and  E^oid  are  the  energy  consumed  to 
load  and  hold  a  registry  value  in  a  shift  register.  Note  that  the  value  of  E^^^  is  defined  based  on  the 
assumption  that  a  Hamming  code  is  exploited  for  the  FEC  which  can  be  implementeded  by  using 
shift  registers  as  well  as  XOR  and  AND  logic  gates  as  described  in  [85].  Therefore,  the  equivalent 
optimization  problem  for  the  FEC  mode,  is  given  by  [PI],  while  replacing  the  constraints  (81)  and 
(84)  by  (93)  and  (95),  respectively,  and  using  constraint  (88)  instead  of  (85)  and  (86),  while  the 
definition  of  Ef^  is  given  by  (99). 

Error  Prevention  Codes  (EPC)  Constraints  Instead  of  correcting  channel  errors  or  just  detecting 
them  and  asking  for  retransmissions  a  posteriori,  the  idea  of  preventing  channel  errors  from  occurring 
in  advanced  or  a  priori,  has  been  recently  proposed  [85].  The  fundamental  idea  is  that,  by  adjusting 
the  coding  weight,  i.e.,  the  average  number  of  logical  ‘T"s  and  “0"s,  the  molecular  absorption 
noise  at  THz-band  frequencies  and  the  multi-user  interference  in  TS-OOK  can  be  minimized.  This 
effectively  results  into  lower  BERs.  Different  EPC  techniques  are  discussed  in  detail  in  [85]. 

In  this  case,  the  packet  error  rate  is  given  by 

=  1  -  (1  -  (100) 

where  n  is  the  number  of  blocks  per  data  packet  payload,  and  k  refers  to  the  block  size.  Moreover, 
the  packet  round-trip  time  in  EPC  is  defined  as  follows: 


rpEPC  _  rpEPC  ,  rpEPC  ,  ^EPC  ( rp  rpEPC  ^ 

^tx  -  ^  tx,data  ^  code  P s,data^^ Prop  +  ^  decode^ 


(101) 


where  is  the  data  transmission  time  and  is  given  by  //r  where  I  is  the  total  length  of  the 

transmitted  data  and  is  equal  to  {Ldata  +  E^^^),  and  is  the  length  of  the  redundant  bits  added  by 
error  prevention  codes.  and  refer  to  latency  caused  by  coding  and  decoding  processes 

of  the  data  respectively  and  are  given  as  follows  [85]: 


jEPC  _  jEPC 
^ code  ^  decode 


=  '^nTcik, 


(102) 


P^data  packet  success  rate  and  can  be  obtained  from  (100).  Finally,  is  given  as 

follows: 


=  l.l{2Tprop  +  Tl^f,J.  (103) 

Eventually,  the  following  equation  describes  the  energy  required  to  accomplish  the  transmission 
process  Ef^  in  EPC: 


pEPC  _  pEPC  ,  pEPC 
^tx  ^tx,data  ^code  ’ 


(104) 


where  refers  to  the  energy  required  to  transmit  the  data  packet,  and  is  given  by  {I  ■  p^  ■ 

and  E^^J^  stands  for  the  consumed  energy  caused  by  computational  process  of  coding  the  data 
in  transmitter,  and  is  given  by  n  ■  {Id  +  k){Eioad  +  Ehoid),  where  Eioad  and  Ehoid  are  the  energy 
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consumed  to  load  and  hold  a  registry  value  in  a  shift  register,  and  Id  is  the  length  of  useful  data  bits 
in  a  transmitted  block.  Note  that  the  value  of  is  defined  based  on  exploiting  logic  gates  and 
parallel-load  shift  registers  to  implement  the  EPC  as  described  in  [85].  Therefore,  the  equivalent 
optimization  problem  for  the  EPC  mode,  is  given  by  [PI],  while  replacing  the  constraints  (81)  and 
(84)  by  (100)  and  (101),  respectively,  and  using  constraint  (88)  instead  of  (85)  and  (86),  while  the 
definition  of  is  given  by  (104). 

7.2.2  Problem  Solution  Approach  and  Algorithm 

Solution  of  The  Optimization  Problem  [PI]  To  solve  the  optimization  problem  [PI],  we  start 
with  the  equality  constraints  (80),  (81),  and  (82).  As  it  is  mentioned  earlier  in  Sec.  7.2.1,  the  BER 
does  not  depend  on  the  optimization  variable  and  the  constraint  (82)  can  be  defined  as  a  given 
parameter  to  our  optimization  problem.  Moreover  the  equalities  (80)  and  (81)  can  be  merged  with 
the  objective  function  of  [PI].  Eor  the  equality  constraint  (83)  which  contains  the  non-smooth 
maximum  function,  we  can  also  plug  it  in  the  objective  function  by  defining  two  auxiliary  efficiency 
functions  miEdata)  and  miEdata)  as  follows: 


^  tj  Ldatai^-lp^) 

^ly^data)  —  , 

r  0^ 

(105) 

X  1  Edatai^  ~ 

mK^data)  —  •  p  , 

r 

(106) 

Hence  we  can  define  our  objective  function  as  follows: 

7]  =  minimiEdata),  miEdata)) ■ 

(107) 

Note  that  Tji{Ldata)  and  r]2{Ldata)  only  depend  on  the  Ldata  as  the  optimization  variable,  and 
all  other  parameters  are  assumed  to  be  given.  Also  the  functions  if/^  ,6^ ,  and  are  all  functions  of 
Edata  which  are  defined  in  Sections  7.2.1  through  7.2.1  for  three  different  error-control  techniques. 
Now  by  merging  the  constraints  (80),  (81),  (82),  and  (83)  with  the  objective  function  of  [PI],  we 
can  define  the  following  equivalent  optimization  problem: 

Optimization  Problem  [P2] 


Given  : 
bind  : 
Minimize  : 
Subject  to  : 


L* 

data 

-7]  =  max{-T]l{Ldata),  -miEdata)) 
Edata  ^  0 


where  ri\{Ldata)  and  miEdata)  are  given  in  (105)  and  (106)  respectively.  Note  that  here  we  are 
using  the  standard  minimization  problem  by  changing  the  sign  of  the  functions  rj,  Tji{Ldata)  and 
miEdata),  and  using  max  function  instead  of  min.  It  can  be  shown  that  [P2]  is  a  Quasi-Convex 
optimization  problem.  More  specifically,  -miEdata)  and  -miEdata)  are  Quasi-Convex  functions 
and  since  “nonnegative  weighted  maximum"  function  preserves  Quasi-Convexity  [86]  therefore  the 
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objective  function  -ri{Lciata)  is  Quasi-Convex.  Moreover,  the  inequality  constraint  (79)  is  a  convex 
set,  hence  [P2]  is  a  Quasi-Convex  optimization  problem.  For  solving  this  optimization  problem  we 
define  the  epigraph  form  of  the  problem  and  use  a  bisection  method  described  in  Sec.  7.2.2. 

Epigraph  Form  and  Bisection  Algorithm  to  Solve  [PI]  To  solve  the  optimization  problem 
stated  in  [P2]  we  can  define  the  epigraph  form  of  the  problem  as  follows: 

Epigraph  form  of  Optimization  Problem  [P2] 


Given  : 
Find  : 
Minimize  : 
Subject  to  : 


L* 

data 

t 

~  ^{J^data)  ~  t  ^  0 
-  Ldata  <  0 


Now,  our  quasi-convex  optimization  problem  can  be  solved  as  a  sequence  of  convex  feasibility 
problems  as  follows: 


Given  : 
Find  : 
Subject  to  : 


r, 

L* 

data 

~  ^ij^data)  ~  t  ^  0 
-  Lclata  <  0 


(108) 


The  above  feasibility  problem  is  convex,  since  all  its  inequality  constraints  are  convex.  Now 
let  us  define  p*  as  the  optimal  value  of  our  optimization  problem  in  [P2].  If  the  problem  (108)  is 
feasible,  then  p*  <  t,  and  if  it  is  not  feasible,  then  p*  >  t.  Therefore  by  using  a  bisection  method  we 
can  solve  the  our  quasi-convex  optimization  problem,  by  solving  a  convex  feasible  problem  in  each 
iteration.  For  the  bisection  method  we  have  to  determine  a  lower  bound  I  and  an  upper  bound  u  for 
the  possible  values  of  p* .  Since  we  are  dealing  with  efficiency  {-p)  as  our  objective  function  so 
we  know  that  it  always  p  6  [0,  Ij.  Hence  we  can  set  the  bounds  to  /  =  -1  and  u  =  0.  However,  it 
can  be  easily  observed  that  min(-p)  >  max{min{pi),  min{p2)}  >  -1,  and  since  pi  and  p2  are  both 
quasi-convex,  differentiable,  and  continuous,  we  can  directly  calculate  the  min  value  of  them  and 
use  it  for  the  lower  bound  as  /  =  max{min(-pi),  min{-p2)}  which  results  in  less  iterations.  The 
steps  of  the  bisection  method  to  solve  [P2]  are  defined  in  Algorithm  1 . 


7.3  Results  and  Discussions 

7.3.1  Numerical  Results  and  Performance  Evaluation 

In  this  section  we  numerically  study  the  performance  of  different  error-control  techniques,  in 
terms  of  maximum  link  throughout  achievable  under  the  energy  harvesting  and  transmission  delay 
constraints. 

In  our  analysis  we  use  the  following  parameters.  We  consider  that  nanosensors  communicate 
by  using  TS-OOK,  hence  every  bit  of  logic  “1"  is  transmitted  as  a  pulse  (which  is  modeled  as  the 
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Algorithm  1  Bisection  Algorithm  for  solving  [P2] 

Set:  I  =  max{min{-rii),min{-ri2)}  using  equations  (105),  (106) 
Set:  u  =  0,  and  6  >  0 
1:  repeat 
2:  t  :=  (I  +  m)/2 

3:  Solve  the  convex  feasibility  problem  (108) 

4:  if  (108)  is  feasible  then 

5:  u  :=  t 

6:  else 

7:  I  :=t 

8:  end  if 

9:  until  u-l  <  e 


derivative  of  a  one-hundred-femtosecond  long  Gaussian  pulse)  with  energy  of  =  1  aJ,  and  the 
logic  “0"  is  transmitted  as  silence  [63].  The  propagation  of  the  pulses  is  modeled  by  utilizing  the 
THz-band  channel  model  in  [23].  Both  TS-OOK  and  this  model  have  been  validated  by  means  of 
extensive  simulations  with  COMSOL  Multi-physics.  Due  to  the  limited  computational  capability  of 
nanomachines,  only  very  simple  error-control  methods  are  utilized.  For  the  ARQ,  we  use  a  16-bit 
CRC  for  error  correction  with  a  2  bytes  long  acknowledgement  packet;  a  Hamming(15,l  1)  code  is 
assumed  for  the  FEC;  and,  for  EPC,  a  16-bit  low-weight  code  with  codeword  size  of  19  bits  is  used 
(EPC  type  II  [85]).  The  resulting  probability  to  transmit  a  pulse  is  0.5  for  ARQ  and  PEC,  and 
0.31  for  EPC. 

We  further  consider  that  the  nanosensors  communicate  with  the  bit  rate  of  r  =  100  Gbps,  and 
the  the  clock  period  to  compute  the  latency  caused  by  the  CRC  or  coding  and  decoding  processes 
is  Tcik  =  1  ps.  We  also  consider  the  energy  required  to  shift,  hold,  and  load  in  a  shift  register  as 
Eshift  =  Ehoid  =  Eload  =  0.1  aJ.  The  communication  distance  range  is  assumed  to  be  J  =  1  -  100 
mm,  and  the  energy  harvesting  rate  ranges  Xharv  =  1  -  400  nJ/s. 

The  link  efficiency  rj  of  the  ARQ  is  shown  in  Pig.  27  as  a  function  of  packet  size  for  a  fixed 
distance  and  energy  harvesting  rate.  in  this  figure  represents  the  link  efficiency  affected  only 
by  the  transmission  time,  i.e.,  the  nano-transmitter  has  enough  energy  to  transmit  and  does  not 
need  time  to  harvest  energy.  In  contrary  shows  the  link  efficiency  which  is  only  affected  by 
energy  harvesting  time  consumption,  i.e.,  the  transmission  time  is  always  less  than  the  time  needed 
to  harvest  energy.  Pinally,  shows  the  trade-off  between  these  two  scenarios  and  shows  the  link 
efficiency  considering  both  constraints.  As  shown  in  Pig.  27,  as  we  increase  the  packet  size,  at  some 
point,  the  energy  harvesting  time  consumption  becomes  dominant  and  restricts  the  efficiency  of  the 
link. 

In  Pig.  28,  the  link  efficiency  rj  of  different  error-control  techniques  is  shown  as  a  function  of 
packet  size  for  a  fixed  distance  and  energy  harvesting  rate.  As  it  can  be  seen  in  this  figure,  EPC 
performs  better  than  the  other  two  techniques  for  very  small  packet  sizes,  and  as  we  increase  the 
packet  size  PEC  outperforms  EPC,  while  ARQ  has  the  lowest  efficiency. 

Pig.  29  and  Pig.  30  depict  the  optimal  link  efficiency  and  the  optimal  packet  size  respectively, 
for  different  error-control  techniques  as  a  function  of  distance  for  a  fixed  energy  harvesting  rate. 
It  can  be  seen  that  for  very  short  distances  EPC  outperforms  the  other  two  error-control  methods, 
while  transmitting  smaller  packets.  However,  PEC  has  a  better  performance  for  longer  distances 
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Fig.  27:  Link  efficiency  for  ARQ,  when  d  =  \  cm. 


Fig.  28:  Link  efficiency  for  different  error-control  techniques,  when  d  =  \  cm  and  Aharv  =  200  nJ/s. 

and  uses  bigger  paekets  for  transmission.  The  ARQ  optimal  paeket  size  is  in  between  the  other  two 
methods,  and  has  the  lowest  link  efficiency. 


Fig.  29:  Optimal  link  efficiency  for  different  error-control  techniques,  when  Aharv  -  200  nJ/s,  and  =  16 
Kbits. 


Finally,  in  Fig.  31,  we  show  the  optimal  link  efficiency  for  different  error-control  methods 
as  a  function  of  energy  harvesting  rate,  for  a  fixed  distance.  This  figure  shows  that  FEC  and  EPC 
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Fig.  30:  Optimal  packet  size  for  different  error-eontrol  teehniques,  when  Aharv  -  200  nJ/s,  and  =  16 
Kbits. 


have  better  performanee  for  low  energy  harvesting  rates,  while  ARQ  outperforms  the  other  two 
teehniques  when  the  nanomaehines  is  eapable  to  harvest  energy  at  higher  rates. 


Fig.  31:  Optimal  link  efheieney  for  different  error-eontrol  teehniques,  when  d  -  \  mm,  and  =  16  Kbits. 
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8  LINK  LAYER  SYNCHRONIZATION  AND  MEDIUM  ACCESS  CONTROL 
8.1  Motivation 

In  parallel  to  the  development  of  THz  deviees  and  physieal  layer  meehanisms,  there  is  a  need 
to  investigate  new  networking  solutions  for  very-high-speed  wireless  data  networks.  Traditional 
MAC  protoeols  need  to  be  revised  in  light  of  the  properties  of  THz -band  eommunieation.  In 
partieular,  the  THz  band  provides  deviees  with  a  very  large  bandwidth  and,  thus,  these  do  not  need 
to  aggressively  eontend  for  the  ehannel.  In  addition,  sueh  very  large  bandwidth  results  in  very  high 
bit-rates  and,  eonsequently,  very  short  transmission  times,  whieh  further  minimize  the  eollision 
probability.  However,  due  to  the  low  transmission  power  of  THz  transeeivers  and  the  high  path-loss 
at  THz -band  frequeneies,  very  high  direetivity  antennas  are  needed  to  establish  wireless  links  beyond 
one  meter.  While  direetional  eommunieation  further  deereases  multi-user  interferenee,  it  requires 
tight  synehronization  between  transmitter  and  reeeiver  to  overeome  the  deafness  problem  [87]. 
Moreover,  the  relatively  long  propagation  delay  when  transmitting  at  Tbps  over  multi-meter-long 
links  results  in  low  ehannel  utilization.  While  these  two  issues  vanish  in  nanoseale  applieations, 
there  are  other  ehallenges  that  affeet  the  link  layer  design  of  THz  nano-deviees,  sueh  as  their  temporal 
energy  fiuetuations  and,  thus,  availability,  introdueed  by  required  energy  harvesting  systems  [88]. 

To  the  best  of  our  knowledge,  there  are  no  existing  protoeols  for  maeroseale  THz-band 
eommunieation  networks.  Available  solutions  for  lower-frequeney  systems  eannot  direetly  be 
utilized  in  this  paradigm,  mainly  beeause  they  do  not  eapture  the  peeuliarities  of  the  THz-band 
ehannel  or  the  eapabilities  of  THz  deviees.  In  terms  of  frequeney,  millimeter-wave  (mm-wave) 
systems  [4],  sueh  as  those  at  60  GHz,  are  the  elosest  existing  teehnology.  In  addition  to  the  IEEE 
802. 1 1  ad  standard  [89] ,  whieh  mainly  adopts  a  very  similar  link  layer  to  that  of  the  entire  IEEE  802. 1 1 
family,  a  few  alternative  MAC  protoeols  for  mm-wave  systems  have  already  been  developed  [90-92]. 
These  protoeols  are  mainly  aimed  at  solving  the  deafness  problem  introdueed  by  direetional  antennas. 
As  in  many  of  the  direetional  MAC  protoeols  for  lower  frequeney  bands  [93,94],  the  existing 
solutions  assume  that  a  node  ean  alternate  between  direetional  and  omnidireetional  antenna  modes. 
However,  at  THz-band  frequeneies,  highly-direetional  antennas  are  simultaneously  needed  both  at 
the  transmitter  and  the  reeeiver  to  sueeessfully  establish  the  link. 

As  part  of  this  effort,  we  have  developed  a  synehronization  and  MAC  protoeol  for  very-high- 
speed  wireless  eommunieation  networks  in  the  THz  band.  The  protoeol  relies  on  a  reeeiver-initiated 
handshake  as  a  way  to  guarantee  synehronization  between  transmitter  and  reeeiver.  In  addition, 
it  ineorporates  a  sliding  window  flow  eontrol  meehanism,  whieh  eombined  with  the  one-way 
handshake,  maximizes  the  ehannel  utilization.  We  eonsider  two  different  applieation  seenarios:  a 
maeroseale  seenario,  in  whieh  nodes  utilize  high-speed  turning  direetional  antennas  to  periodieally 
sweep  the  spaee,  and  a  nanoseale  seenario,  in  whieh  nodes  make  use  of  a  piezoeleetrie  nano-generator 
to  harvest  energy.  Eor  eaeh  seenario,  we  eonsider  a  different  physieal  layer,  namely,  a  traditional 
earrier-based  modulation  for  the  maero-seenario  and  a  femtoseeond-long  pulse -based  modulation 
with  user  interleaving  for  the  nano-seenario.  We  analytieally  investigate  the  performanee  of  the 
proposed  MAC  protoeol  in  terms  of  delay,  throughput  and  sueeessful  paeket  transmission  probability, 
and  eompare  it  to  that  of  a  modified  Carrier  Sense  Multiple  Aeeess  with  Collision  Avoidanee 
(CSMA/CA)  with  and  without  handshake.  In  addition,  we  have  implemented  in  ns-3  the  proposed 
protoeol  and  the  neeessary  THz  models  (ehannel,  earrier-based  and  pulse -based  physieal  layers, 
turning  antenna  and  harvesting  system),  and  provide  extensive  simulation  results  to  validate  the 
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performance  of  our  solution. 

Next,  in  8.2.1,  we  describe  the  system  model  considered  throughout  the  paper  and  derive 
formulations  for  the  bit  and  frame  error  rate,  the  probability  of  collision  in  the  macroscale  and  the 
nanoscale  scenarios,  and  the  energy  harvesting  system  for  nano-devices.  In  Sec.  8.2.2,  we  describe 
the  proposed  protocol  and  analytically  investigate  its  performance  for  the  two  scenarios.  We  provide 
simulation  and  numerical  results  in  Sec.  8.3. 

8.2  Methods,  Assumptions  and  Procedures 
8.2.1  THz-band  Communication  System  Model 

In  this  section,  we  summari/e  the  main  peculiarities  of  THz-band  communication  networks,  both 
for  the  macroscale  and  the  nanoscale  scenarios.  In  particular,  first,  we  formulate  the  Signal-to-Noise 
Ratio  (SNR),  which  is  needed  for  the  computation  of  Bit  Error  Rate  (BER)  and  Erame  Error  Rate 
(EER),  starting  from  an  accurate  THz-band  channel  model.  Then,  we  formulate  the  collision 
probability  in  the  macroscale  scenario,  for  which  we  introduce  and  analyze  the  impact  of  utilizing 
high-speed  turning  THz  directional  antennas,  as  well  as  in  the  nanoscale  scenario,  for  which  we 
revise  the  concept  of  interleaved  pulse-based  transmissions  [95].  Einally,  we  introduce  the  energy 
harvesting  model  utilized  in  nanonetworks. 

Signal- to- noise  Ratio,  Bit  Error  Rate  and  Frame  Error  Rate  The  propagation  of  electromag¬ 
netic  waves  at  THz-band  frequencies  is  mainly  affected  by  molecular  absorption,  which  results 
in  both  molecular  absorption  loss  and  molecular  absorption  noise.  In  particular,  based  on  the 
THz-band  channel  model  introduced  in  [23],  the  signal  power  at  a  distance  d  from  the  transmitter, 
Pr  is  given  by 

Pr  (d)  =  f  St  if)  \H,  if,  d)|2  \Hr  (/)|2  df,  (109) 

Jb 

where  St  is  the  single-sided  power  spectral  density  (p.s.d)  of  the  transmitted  signal,  B  stands  for  its 
bandwidth  and  /  refers  to  frequency.  He  refers  to  the  THz-band  channel  frequency  response,  which 
is  given  by 

where  c  refers  to  the  speed  of  light  and  kabs  is  the  molecular  absorption  coefficient  of  the  medium. 
This  parameter  depends  on  the  molecular  composition  of  the  transmission  medium,  i.e.,  the  type 
and  concentration  of  molecules  found  in  the  channel  and  is  computed  as  in  [23].  Hr  in  (109)  refers 
to  the  receiver  frequency  response,  which  we  consider  an  ideal  low-pass  filter  with  bandwidth  B,  for 
the  time  being. 

Similarly,  the  molecular  absorption  noise  power  Nr  at  a  distance  d  from  the  transmitter,  which 
can  be  modeled  as  additive,  Gaussian,  colored  and  correlated  to  the  transmitted  signal  [63],  is  given 
by 

Nr  (d)  =  f  (S^B  (/)  +  S^l  (/,  d))  \Hr  (/)|2  df,  (1 1 1) 

Jb 

where  it  is  taken  into  account  that  the  total  molecular  absorption  noise  is  contributed  by  the 
background  atmospheric  noise  p.s.d.,  Sj^b  and  the  self-induced  noise  p.s.d.,  S^/,  and  are  computed 
as  described  in  [63]. 
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The  SNR  at  a  distance  d  from  the  transmitter  can  be  then  obtained  dividing  (109)  by  (1 1 1). 
For  a  given  SNR,  the  BER,  Pt,,  only  depends  on  the  modulation  technique.  In  our  analysis,  we 
consider  two  different  physical  layers: 

•  In  the  macroscale  scenario,  i.e.,  for  distances  above  a  few  meters,  we  consider  that  nodes 
transmit  a  conventional  m-PSK  modulated  signal  over  a  100-GHz-wide  transmission  window 
at  1.05  THz.  This  is  done  to  overcome  the  very  high  molecular  absorption  loss  over  long 
distances  [23]. 

•  In  the  nanoscale  scenario,  i.e.,  for  distances  much  below  one  meter,  we  consider  that  nodes 
transmit  by  using  TS-OOK,  a  modulation  scheme  based  on  the  transmission  of  one-hundred- 
femtosecond-long  pulses  by  following  an  on-off  keying  modulation  spread  in  time  [63].  The 
p.s.d.  of  such  pulses  is  mainly  contained  between  0.9  and  4  THz. 

Finally,  to  compute  the  frame  error  rate  Pp,  we  consider  that  a  hybrid  error  control  technique 
based  on  the  combination  of  low- weight  channel  codes  with  ARQ  is  utilized  [67] . 

Collision  Probability  The  collision  probability  depends  on  the  application. 

A.  Macroscale  Scenario  In  this  case,  nodes  require  high-gain  directional  antennas  in  transmission 

and  in  reception  to  successfully  establish  a  link.  As  we  already  discussed,  this  results  into  low 
multi-user  interference,  at  the  cost  of  high  synchronization  requirements.  To  model  the  multi-user 
interference,  we  need  to  take  into  account  both  the  spatial  distribution  of  the  nodes  as  well  as 
their  temporal  activity.  In  our  model,  we  consider  that  nodes  are  randomly  distributed  in  space 
by  following  a  spatial  Poisson  process  with  rate  Aa-  The  area  of  influence  of  an  individual  node 
is  given  by  A(0)  =  =  p,  where  9  is  the  antenna  beam- width  in  radians  and  stands  for  the 

maximum  transmission  distance.  Then,  the  probability  of  Ending  i  nodes  in  A  is  given  by 

P  [i  6  A(0)]  =  g-^AA(e)^  ^2) 

i\ 

If  each  node  i  in  A  generates  new  frames  at  a  rate  given  by  1  jaFT,  where  a  is  a  constant  and  FT 
stands  for  frame  time,  the  aggregated  traffic  generated  by  i  nodes  is  simply  Aj  =  ijaFT.  Thus,  the 
probability  that  j  out  of  i  nodes  are  active  during  a  vulnerable  time  of  2FT  is  given  by: 

P  [j  e  2FT]  =  ^ ^  ^3^ 

j! 

Finally,  the  collision  probability  in  the  macroscale  scenario  is: 

00 

P,  =  J]p[ieA  (91)]  (1  -  P  [0  6  2FT]) .  (1 14) 

i=l 

B.  Nanoscale  Scenario  In  this  case,  miniature  nodes  do  not  require  directional  antennas,  but 
transmit  omnidirectionally.  Therefore,  the  probability  of  finding  i  nodes  in  the  area  of  influence  of  a 
node  is  given  by  (1 12)  for  0  =  2n.  This  results  into  a  much  larger  number  of  potentially  interfering 
nodes.  Interestingly,  however,  TS-OOK  supports  the  simultaneous  transmission  and  reception  of 
time -interleaved  frames.  In  this  scheme,  a  logical  “0"  is  transmitted  as  silence,  whereas  a  logical 
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“1"  is  transmitted  with  a  pulse.  The  time  between  symbols,  pulses  and  silenees,  Ts  is  mueh  longer 
than  the  pulse  duration  Tp,  /3  =  Ts/Tp  »  1. 

By  eonsidering  that  eaeh  node  i  in  A  generates  new  frames  at  a  rate  1  jaFT,  the  rate  at  whieh 
new  pulses  are  generated  is  given  hy  Xp  =  ip\  jaTs,  where  p\  refers  to  the  probability  of  transmitting 
a  pulse,  and  is  related  to  the  eoding  weight  [67].  Then,  the  probability  of  j  out  of  i  nodes  generate 
pulses  within  the  vulnerable  time  ITp  is: 


P  [;  e  ITp]  =  i^f^ELe-^p2Tp  _  (115^ 

Finally,  if  there  are  n  symbols  in  one  frame,  the  eollision  probability  of  one  frame  in  the  nanoseale 
seenario  is: 

00 

Pc  =  J]P[ie  A(2n)]  (l  -F  [O  6  2Tp]"^  .  (116) 

i=l 

Nano  Energy  Model  In  our  analysis,  we  are  interested  in  eomputing  the  frame  error  probability 
eaused  by  insuffieient  energy  at  the  transmitter  or  at  the  reeeiver.  For  this,  we  utilize  the  same 
methodology  as  in  [36].  Due  to  spaee  limitations,  we  only  qualitatively  deseribe  our  methodology, 
but  refer  the  reader  to  the  original  paper  for  further  details.  Mainly,  the  probability  mass  funetion 
(p.m.f.)  of  the  energy  stored  at  the  nano-node  battery  after  reaehing  a  steady  state  depends  both  on 
the  rate  at  whieh  energy  is  harvested  and  the  rate  at  whieh  energy  is  eonsumed.  The  latter  depends  on 
the  the  new  frame  generation  rate  as  well  as  on  the  expeeted  number  of  retransmissions.  Ultimately, 
different  MAC  protoeols  result  into  different  number  of  transmissions  and,  thus,  into  different  p.m.f. 
for  the  energy  at  the  nano-battery.  With  the  p.m.f.  of  the  energy,  we  ean  ealeulate  the  probability  of 
having  enough  energy  at  eaeh  node.  We  implement  the  harvesting  model  in  ns-3  and  utilize  the 
eolleeted  data  to  estimate  the  battery  p.m.f.  used  in  our  numerieal  results. 

8.2.2  Receiver-initiated  Synchronization  and  Medium  Access  Control  Protocol 

In  this  section,  we  first  describe  the  proposed  protocol  for  THz-band  communication.  Then,  we 
analytically  investigate  its  performance  in  terms  of  successful  packet  delivery  probability,  packet 
delay  and  throughput. 

Protocol  Overview  In  light  of  the  system  model  described  in  Sec.  8.2.1,  it  is  clear  that  an  initial 
handshake  is  needed  in  THz-band  communication  networks  to  guarantee  link-layer  synchronization 
between  the  transmitter  and  the  receiver.  The  objective  of  such  handshake  is  to  prevent  unnecessary 
data  transmissions  when  the  receiver  is  not  available,  whether  because  it  is  not  facing  the  transmitter 
(macroscale  scenario)  or  because  it  does  not  have  enough  energy  to  handle  a  new  transmission 
(nanoseale  scenario).  The  fundamental  idea  behind  the  proposed  protocol  is  to  reduce  the  overhead 
introduced  by  such  handshaking  process  by  having  nodes  announce  their  availability  to  receive 
data.  In  other  words,  the  traditional  two-way  handshake  is  reduced  to  a  one-way  handshake  process. 
Receiver-initiated  MAC  protocols  have  been  successfully  utilized  in  other  scenarios  [93,94,96], 
but  the  existing  solutions  cannot  directly  be  utilized  in  our  scenario  because  of  the  aforementioned 
peculiarities  of  THz-band  communication  networks. 

Besides  the  one-way  handshake,  the  proposed  protocol  also  makes  use  of  sliding  flow  control 
window  at  the  link  layer  to  maximize  the  channel  utilization.  In  particular,  both  the  delay  introduced 
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while  waiting  for  the  receiver  availability  and  the  relatively  long  propagation  delay  when  transmitting 
at  multi-Gbps  or  Tbps  over  multi-meter  distances  (macroscale  scenario)  result  in  a  relative  low 
channel  utilization.  To  overcome  such  problem,  the  receiver  can  specify  the  time  that  it  will  remain 
facing  in  the  current  direction  (macroscale  scenario)  or  the  amount  of  data  packets  willing  to  accept 
with  its  currently  available  energy  (nanoscale  scenario). 

The  basics  of  the  protocol  are  summarized  next.  A  node  can  be  found  in  transmitting  mode  or 
in  receiving  mode: 

•  A  node  in  transmitting  mode  (TM),  i.e.,  with  data  to  transmit,  checks  whether  a  current 
Clear-To-Send  (CTS)  frame  from  the  intended  receiver  has  been  recently  received.  We 
consider  that  a  CTS  has  an  expiration  time,  which  is  a  parameter  value  in  our  system.  If  not, 
the  node  listens  to  the  channel  until  the  reception  of  a  new  CTS  frame.  In  the  macroscale 
scenario,  we  consider  that  the  node  in  TM  points  its  directional  antenna  to  the  receiver.  For 
this,  we  consider  that,  at  the  link  layer,  the  node  in  TM  knows  the  position  of  the  receiver. 
This  information  is  provided  by  the  network  layer  after  a  discovery  process,  which  remains 
as  future  work.  In  the  nanoscale  scenario,  we  consider  that  nodes  utilize  omnidirectional 
antennas. 

•  A  node  in  receiving  mode  (RM),  i.e.,  with  sufhcient  resources  (e.g.,  energy  or  even  memory) 
to  handle  a  new  incoming  transmission,  broadcasts  its  status  by  means  of  a  CTS  frame.  In  the 
macroscale  scenario,  the  node  in  RM  uses  a  dynamically  turning  narrow-beam  to  broadcast 
CTS  frames  while  sweeping  its  entire  surrounding  space.  Such  electrically-controlled  high¬ 
speed  turning  directional  antennas  can  be  implemented  for  example  by  means  of  the  proposed 
graphene-based  plasmonic  nano-antenna  array  (Sec.  4).  The  node  in  RM  mode  cannot  know 
in  advance  who  is  willing  to  transmit,  which  is  why  it  needs  to  sweep  the  entire  space.  The 
turning  speed  of  beam  is  a  parameter  to  be  optimized  in  our  analysis.  In  the  nanoscale  scenario, 
this  information  is  omnidirectionally  transmitted.  The  CTS  frame  also  contains  information 
on  the  receiver’s  current  sliding  window  size. 

Upon  the  reception  of  a  CTS  frame,  a  node  in  TM  mode  checks  whether  it  has  data  for  such 
receiver  and  the  necessary  resources.  If  so,  it  proceeds  with  the  DATA  frame  transmission  by  taking 
into  account  also  the  receiver  sliding  window  size.  If  the  transmission  is  successful,  the  node  in  RM 
sends  a  positive  acknowledgement  (ACK)  frame.  Otherwise,  after  a  time-out,  the  node  in  TM  will 
set  a  random  back-off  time,  which  depends  on  the  number  of  transmission  attempts,  and  repeats  the 
entire  process  when  done.  After  successfully  receiving  a  packet,  i.e.,  successfully  transmitting  a 
CTS,  a  DATA  and  an  ACK  frames,  the  node  in  RM  can  decide  to  keep  turning,  continue  to  collect 
more  packets,  or  switch  to  TM. 

A  couple  of  comments  regarding  fairness  need  to  be  made.  First,  as  in  any  triggered  reaction 
protocol,  nodes  in  TM  wait  a  random  time  after  receiving  a  CTS  frame  and  before  sending  the  DATA 
frame.  In  the  macroscale  scenario,  carrier-sense  is  performed  during  that  time.  In  the  nanoscale 
scenario,  in  which  a  pulsed-based  physical  layer  is  used,  there  is  no  carrier  to  sense,  but  the  chances 
of  having  a  collision  are  very  low  (Sec.  8.2.1).  Second,  only  for  the  macroscale  scenario,  it  is 
relevant  to  note  that  the  node  in  RM  cannot  simply  stop  indefinitely  at  a  node  in  TM,  but  it  needs  to 
continue  “turning".  Therefore,  we  need  to  guarantee  that,  within  the  small  span  of  time  that  the  node 
in  RM  is  facing  the  node  in  TM,  the  DATA  frame  can  be  successfully  transmitted.  This  is  possible 
because  nodes  are  transmitting  at  multi-Gbps  and  even  Tbps  and,  thus,  only  several  nanoseconds  are 
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usually  needed,  whieh  is  mueh  lower  than  in  existing  wireless  eommunieation  systems  and,  thus, 
ean  be  leveraged  to  effeetively  never  stop  turning. 

Performance  Analysis  To  analyze  the  performanee  of  the  proposed  protoeol,  we  analytieally 
investigate  the  sueeessful  paeket  delivery  probability,  the  paeket  delay  and  the  throughput.  First, 
we  foeus  on  the  maeroseale  seenario  and  develop  a  mathematieal  framework  in  detail.  Then,  we 
summarize  the  differenees  needed  to  eapture  the  nanoseale  seenario  peeuliarities. 

A.  Macroscale  Scenario  In  this  case,  the  main  factor  affecting  the  performance  of  the  protocol  is 
the  antenna  turning  speed  co,  given  in  circles-per-second.  In  particular,  we  consider  that  the  antenna 
shifts  its  direction  in  discrete  steps  and,  thus,  provides  coverage  to  different  sectors  in  different  times. 
We  define  the  sector  time  or  time  during  which  the  antenna  beam  is  pointing  to  a  certain  direction  as 
Tsector  =  0l{2nco).  This  limits  the  maximum  number  of  retransmissions  that  a  node  can  complete  in 
the  current  round  s.  This  affects  the  overall  packet  delay  and  throughput,  as  the  node  will  have  to 
wait  for  an  entire  cycle  before  being  able  to  continue  its  ongoing  transmission  in  round  5  +  1 .  More 
specifically,  the  maximum  number  of  retransmission  rjmax  that  a  node  can  complete  in  round  s  can 
be  calculated  as  follows 


r]max  [5]  =  min 


t  sector 


Tcts  -  T, 


prop 


Tt/o  +  Thjc 


(117) 


where  =  ITp^op  +  +  Tdata  +  Tack  is  the  time-out  time,  Tcts,  Tdata  and  Tack  refer  to  the 

CTS,  DATA,  and  ACK  frames  transmission  time,  respectively,  and  Ti^jo  is  a  random  exponential 
back-off  time.  In  our  analysis,  we  do  not  ignore  the  impact  of  the  propagation  delay  Tprop,  as  it 
is  comparable  or  even  larger  than  the  frame  transmission  time  in  the  macroscale  scenario,  k  is  a 
parameter  value  that  specifies  the  maximum  number  of  retransmissions  available  in  the  current 
round.  In  particular,  k  [1]  is  set  to  a  default  value  ko.  For  example,  if  ko  =  5,  the  total  amount 
of  retransmissions  for  a  specific  packet  is  equal  to  5.  This  can  be  “consumed"  within  one  round 
if  Tsector  is  Very  long,  i.e.,  the  antenna  turns  at  a  slow  speed,  or  might  be  spread  across  rounds, 
otherwise. 

Then,  the  probability  to  succeed  with  exactly  i  retransmissions  within  the  same  round  is  given 

=  PcTS  (1  -  PdataPackY  ^  PdataPack,  (118) 

where  Fcrs  =  =  (1  -  Fp)  is  the  probability  of  successfully  receiving  a  CTS  frame;  =  PcPp 

is  the  probability  of  successfully  receiving  a  DATA  frame;  Pack  =  Pp  is  the  probability  of 
successfully  receiving  an  ACK  frame.  We  consider  that  the  main  reason  for  not  properly  receiving 
the  CTS  is  the  presence  of  bit  errors,  rather  than  the  collision  with  other  CTS.  In  the  macroscale 
scenario,  this  is  generally  true,  as  it  is  very  unlikely  to  have  two  or  more  receivers  exactly  pointing 
towards  the  same  transmitter  with  their  directional  beams.  Similarly,  in  the  nanoseale  scenario, 
given  the  energy  constraints  of  nano-nodes,  it  is  not  likely  to  have  many  nearby  receivers  announcing 
their  availability  at  the  same  time.  When  it  comes  to  the  DATA  frame,  failures  might  occur  because 
of  both  bit  errors  as  well  as  collisions.  Even  if  we  introduce  a  random  initial  delay  between  the  CTS 
reception  and  the  DATA  transmission,  collisions  can  occur.  Finally,  ACK  frames  might  also  suffer 
from  bit  errors. 

From  this,  we  can  easily  write  the  probability  to  successfully  transmit  the  packet  in 
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round  as  well  as  the  expected  number  of  retransmissions  rj  in  that  round  as 


p::iM=  E  r,is]=  iPiiSf-  dw) 

i=l  !=1 

If  successful,  the  average  successful  packet  delay  introduced  by  the  current  round  5  can  be  calculated 
as 

Tsucc  [5]  =  (V  [5]  -  1)  {Tt/o  +  Tb/o)  +  T’lco  (120) 

where  T}^cc  -  '^Tprop  +  Tproc  +  To  at  a  +  Tack  is  the  delay  when  successfully  transmitting  the  packet 
in  one  attempt. 

If  the  node  is  not  successful  in  the  current  round,  but  the  maximum  number  of  retransmissions 
across  rounds  ko  has  not  been  yet  achieved,  the  node  waits  for  a  new  CTS  frame  (after  one  antenna 
cycle).  Now  k  [5  +  1]  =  k  [5]  -  77™“*  [5].  The  maximum  number  of  rounds  Tmax  is  given  by 

T max  =  min  {k  [5]  =  0}  .  (121) 


Then,  the  probability  to  successfully  transmit  a  packet  in  the  j-th  round  is  given  by 


pj-rnd 
^  succ 


\-p 


nid 


prnd 
^  succ 


[j]. 


(122) 


where  Psucc  1®  given  in  (1 19). 

From  this,  the  total  packet  successful  delivery  probability  and  the  average  number  of  rounds 
needed  to  do  so  are  given  by 


r  r 

^  max  ^  max 

Psucc  =  £  PiZc"-,  r  =  £  jPiZf.  (123) 

7=1  7=1 

The  discard  probability  immediately  follows  as  Pdis  =  Psucc-  The  average  packet  delay  can  similarly 
be  obtained  as: 


^  packet 


^max  j  ^ 

=En( 


1  -  [«])  m 


i=\  u=\ 

((j  ~  T)Tcycle  "I"  Tsucc  [7])  Tujaitt 


(124) 


where  Tcyde  is  the  time  needed  for  the  antenna  to  complete  one  entire  circle.  T^ait  refers  to  the 
average  time  that  the  transmitter  will  have  to  wait  for  the  receiver’s  CTS  in  the  first  round,  and  is 
computed  as  follows 

^sectors 

Twait  ~  E  p/~^Pf(i-l)T,ector,  (125) 

i=l 

where  n sectors  =  2;7r/0  is  the  number  of  sectors  and  Pf  =  6  jn  the  transmitter  and  receiver  facing 
probability.  Finally,  we  can  obtain  the  throughput  as  5  =  n  packet  IT  packet,  where  ripacket  is  the  number 
of  bits  per  packet. 
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B.  Nanoscale  Scenario  In  this  case,  the  main  problem  affecting  the  performance  of  the  protocol 
is  posed  by  the  fluctuations  in  the  available  energy  in  each  nano-device,  which  were  discussed  in 
Sec.  8.2.1.  In  our  analysis,  we  consider  that  the  battery  in  each  nano-node  can  hold  up  to  emax 
energy  frames.  An  energy  frame  is  the  energy  in  Joules  consumed  in  the  reception  of  a  control  frame 
(CTS  or  ACK),  also  denoted  as  6'’^  ,  ,.  Similarly,  ,  ,  denotes  the  energy  frames  required  to 
transmit  a  control  frame,  and  and  are  the  number  of  energy  frames  needed  to  transmit 
and  receive  a  DATA  frame,  respectively.  The  rate  at  which  energy  frames  are  harvested  is  denoted 
by  Aharv  Based  on  the  probability  distribution  of  the  battery  energy  status,  we  can  calculate  the 
insufficient  energy  probability  and  its  impact  on  the  system. 

The  model  now  is  fundamentally  the  same,  with  only  few  modifications.  In  this  case,  the 
reception  of  a  CTS  frame  is  not  governed  by  the  antenna  turning  speed  co  and  the  sector  time  T sector, 
but  by  the  time  needed  by  the  receiver  to  harvest  enough  energy  to  operate  and  the  announced  CTS 
lifetime.  Upon  the  reception  of  a  CTS  frame,  the  transmitter  node  checks  whether  it  has  enough 
energy  to  successfully  transmit  one  or  more  packets.  The  new  maximum  number  of  retransmissions 
ijniax  in  round  5  is  given  by 


rimax  [■^]  =  niin 


CTSlife 

Tfio  +  Thjo  J 


,k[s]. 


'-max 

^  P  [level  =  /] 


tx 


I  -  e 


L  ^retry  \ 


+  1 


(126) 


where  CTSufe  refers  to  the  CTS  frame  lifetime,  set  by  the  receiver  according  to  its  energy,  Tt/o 
and  are  the  time-out  and  back-off  times,  respectively,  defined  similarly  as  for  the  macroscale 
scenario,  ^’ack  energy  required  to  complete  a  packet  transmission  on  the 

transmitter  side,  level  is  the  number  of  energy  frame  units  contained  in  the  battery,  and  eretry  is  the 
energy  spent  in  a  retransmission.  As  before,  k  is  a  parameter  value  that  specifies  the  maximum 
number  of  retransmissions  still  available. 

In  general,  there  is  no  guarantee  that  the  transmitter  and  the  receiver  will  have  enough  energy 
to  successfully  complete  the  packet  transaction  in  one  round,  i.e.,  during  the  duration  of  the  current 
CTS  lifetime.  As  a  result,  both  nodes  might  have  to  wait  to  harvest  the  required  energy.  The  average 
waiting  times  for  the  transmitter  and  the  receiver  nodes  are  given  by 


rptx/rx 
^  wait 


1 


^ha 


^txjrx 

^  P  [level  =  /] 
!=0 


(127) 


where  -l-  is  the  required  energy  to  start  transmission  from  receiver  side, 

and  Aharv  is  the  energy  harvesting  rate  in  frames  per  second.  Thus,  the  average  waiting  time  for  the 
packet  transmission  can  be  calculated  as: 


wait  ^  wait^  wait^  wait  ^  wait^  wait^  wait  ^  wait^  wait  ^  wait' ^  wait  J  ^ 


(128) 


where 

^tx  ^rx 

Pwait  =  =  '■]  ;  Pwlit  =  J]P  =  '■]  ’  (129) 

i=0  i=0 

are  the  waiting  probability  for  the  transmitter  and  the  receiver,  respectively.  The  packet  delay  Tdeiay 
can  be  now  obtained  by  utilizing  (128)  in  (124)  instead  of  (125).  The  rest  of  the  analytical  model 
remains  the  same. 
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8.3  Results  and  Discussions 


In  this  section,  we  investigate  the  performance  of  the  proposed  protocol  and  compare  it  to  that 
of  CSMA/CA  with  and  without  RTS/CTS  (2-way  and  0-way  handshake).  The  numerical  results 
obtained  by  solving  the  analytical  models  developed  in  Sec.  8.2.2  are  validated  by  simulation 
obtained  with  ns-3.  For  this,  we  have  implemented  the  frequency-selective  THz-band  channel,  the 
two  THz  physical  layers  (carrier-based  and  pulse-based  with  interleaving),  the  high-speed  turning 
antenna,  the  energy  harvesting  unit,  our  proposed  protocol  and  we  also  tailored  and  tuned  CSMA/CA 
to  work  with  the  THz  models.  For  completeness,  we  have  also  developed  the  analytical  models  for 
CSMA/CA  with  and  without  RTS/CTS,  but  these  are  not  included  due  to  space  constraints. 

8.3.1  Macroscale  Scenario 

In  this  case,  we  utilize  the  following  parameter  values.  The  network  is  composed  by  18  nodes  over  a 
circular  area  with  radius  I  =  10  m.  The  transmission  power  of  each  node  is  limited  to  10  /iW.  Nodes 
are  equipped  with  turning  directional  antennas  with  directivity  D  =  35  dB,  which  corresponds  to  an 
antenna  beam- width  of  approximately  32°.  The  center  frequency  fc  and  bandwidth  B  are  1.05  THz 
and  100  GHz,  respectively,  which  corresponds  to  the  hrst  absorption-dehned  transmission  window 
above  1  THz  (any  other  transmission  window  can  be  selected).  The  antenna  turning  speed  is  a 
parameter  in  our  analysis.  The  back-off  time  4/^  and  the  sector  time  T sector,  which  is  utilized  to  also 
set  the  receiver  sliding  window,  is  adjusted  according  to  the  turning  antenna  speed.  The  maximum 
number  of  retransmissions  ko  is  set  to  5. 

In  Fig.  32,  the  packet  discard  probability  is  shown  as  a  function  of  the  antenna  turning  speed. 
As  expected,  with  our  proposed  receiver-initiated  or  1-way  handshake  protocol,  the  probability  of 
discarding  a  packet  is  virtually  zero  and  signihcantly  much  lower  than  0-way  and  2- way  protocols. 
The  main  reason  for  this  is  that  no  retransmission  attempts  will  be  “wasted"  when  the  receiver  is 
not  facing  the  transmitter,  i.e.,  unless  the  transmitter  has  recently  received  a  CTS  frame  from  the 
intended  receiver.  This  is  not  the  case  for  the  two  traditional  protocols,  which  might  discard  up  to 
30%  of  their  generated  packets  as  the  facing  probability  given  in  (125)  is  very  low  for  such  very 
high  directivity  antennas. 

However,  the  cost  of  a  lower  discard  probability  is  reflected  in  the  achievable  throughput, 
which  is  plotted  in  Fig.  33  as  a  function  of  the  antenna  turning  speed.  For  low  turning  antenna 
speeds,  the  throughput  achieved  by  0-way  or  2-way  protocols  is  higher  than  that  of  the  proposed 
protocol,  mainly  because,  successful  packets  (less  than  80%  of  the  total)  have  been  transmitted 
when  the  receiver  has  been  found  facing  at  the  transmitter,  without  the  need  to  wait  for  it.  As  the 
antenna  turning  speed  increases,  the  throughput  for  the  proposed  protocol  increases  and  ultimately 
meets  that  of  the  other  two  protocols  with  the  advantage  of  having  no  packets  dropped.  The  values 
considered  for  the  turning  speed  are  conservative  in  light  of  the  theoretically  achievable  with  the 
aforementioned  structures. 

8.3.2  Nanoscale  Scenario 

In  this  case,  we  utilize  the  following  parameter  values.  A  circular  area  with  radius  I  =  0.01  m 
with  varying  nano-node  densities  is  considered.  Nano-nodes  communicate  by  utilizing  TS-OOK 
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Fig.  32:  Discard  probability  as  a  function  of  the  antenna  turning  speed  o). 


Fig.  33:  Throughput  as  a  function  of  the  antenna  turning  speed  oj. 

with  pulse  energy  Ep=\  attoJoule  (aJ),  pulse  length  rp=100  fs,  and  spreading  faetor  TsjTp  =  100. 
Nodes  utilize  omnidireetional  antennas  with  no  direetivity  gain.  For  the  energy  model,  we  define 
^l^harv  =  8  /ts  per  energy  frame,  ^data=^^  energy  units.  The 

eapaeity  of  the  battery  is  emax=1000  energy  frames.  These  values  have  been  obtained  by  taking  into 
aeeount  the  energy  to  transmit  and  reeeive  a  pulse,  the  eoding  weight,  and  the  harvesting  model 
given  in  [36]. 

In  Fig.  34,  the  paeket  diseard  probability  is  shown  as  a  funetion  of  the  node  density.  Similarly 
as  for  the  maeroseale  seenario,  the  number  of  dropped  paekets  is  signifieantly  lower  than  with  the 
other  protoeols.  It  is  relevant  to  note  that  in  this  ease,  the  0-way  handshake  protoeol  has  also  a 
very  low  diseard  probability.  The  reason  for  this  is  beeause  the  time  needed  to  harvest  energy  to 
transmit  a  paeket  is  mueh  larger  than  the  time  needed  to  harvest  energy  to  reeeive.  As  a  result, 
when  the  transmitter  has  enough  energy  to  transmit,  it  is  likely  that  the  reeeiver  has  it  too.  However, 
for  every  transmission  missed,  the  time  before  the  next  retransmission  is  very  high.  This  ean  be 
seen  in  Fig.  35,  where  the  throughput  is  shown  as  a  funetion  of  the  node  density.  Note  that  in 
this  ease  the  throughput  for  our  proposed  protoeol  is  signifieantly  higher  than  for  the  other  two 
protoeols  (note  the  logarithmie  seale).  The  reason  is  again  related  to  the  very  long  time  needed  for  a 
node  to  have  enough  energy  to  transmit  after  a  failed  transmission.  This  emphasizes  the  need  for 
transmitter-reeeiver  synehronization,  whieh  is  the  main  reason  behind  the  one-way  handshake. 

8.3.3  Discussion 

We  have  presented  a  link-layer  synehronization  and  MAC  protoeol  for  ultra-high- speed  wireless 
eommunieation  networks  in  the  THz  band.  The  protoeol  relies  on  a  reeeiver-initiated  handshake  as 
well  as  a  sliding  window  flow  eontrol  meehanism  to  guarantee  synehronization  between  transmitter 
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Fig.  34:  Packet  discard  probability  as  a  function  of  the  node  density. 


Fig.  35 :  Throughput  as  a  function  of  the  node  density. 


and  receiver,  maximize  the  channel  utilization  and  minimize  the  packet  discard  probability.  The 
performance  of  the  proposed  protocol  is  analytically  investigated,  compared  to  that  of  a  modified 
CSMA/CA  with  and  without  RTS/CTS,  and  validated  though  extensive  simulations  with  ns-2>.  The 
results  show  that  the  proposed  protocol  can  maximize  the  successful  packet  delivery  probability  and 
enhance  the  achievable  throughput  in  THz-band  communication  networks. 
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9  CONCLUSIONS 


The  potential  impaets  of  THz  eommunieation,  onee  realized,  would  be  seeond  to  none.  THz 
teehnology  has  been  reeently  identified  by  DARPA  as  one  of  the  four  major  research  areas  that 
could  eventually  have  an  impact  on  our  society  larger  than  that  of  the  Internet  itself  [97].  The  THz 
band  opens  the  door  to  a  plethora  of  applieations  [8]  in  very  diverse  domains,  ranging  from  Terabit 
Wireless  Personal  and  Loeal  Area  Networks  to  wireless  nanosensor  networks  [29]  or  the  Internet 
of  Nano-Things  [34].  In  this  eontext,  the  development  of  a  new  eommunieation  and  networking 
teehnology  to  support  networks  with  “billions  of  eonneeted  nanosystems”  has  been  reeently  identified 
as  one  of  the  four  essential  components  of  the  next  IT  revolution  by  the  Semiconductor  Research 
Consortium  and  NSF  [98]. 

As  part  of  this  projeet,  we  have  made  novel  eontributions  aeross  different  aspeets  required  to 
enable  praetieal  THz  eommunieation  systems.  For  a  period  of  20  months  and  with  a  total  budget  of 
$240,000,  we  have  developed  both  novel  deviee  eoneepts  as  well  as  physieal  and  link  layer  solutions, 
whieh  resulted  in  a  total  of  eight  eonferenee  paper  publieations  [24, 26, 28, 30, 31, 33, 37, 38],  three 
journal  publieations  [27,32,39],  four  additional  journals  submitted  for  publieation  or  in  preparation, 
and  potentially  one  patent  applieation. 

As  part  of  this  projeet,  four  graduate  students  have  been  partially  supported.  The  developed 
work  has  been  presented  at  eonferenees  by  the  Pis  as  well  as  their  students,  who  have  been  trained  to 
beeome  experts  in  this  quiekly  evolving  field.  Through  their  partieipation  in  eonferenees,  they  have 
also  been  exposed  to  other  researeh  groups  and  researeh  areas,  and  have  now  a  better  understanding 
of  the  broad  wireless  eommunieation  eommunity. 

As  with  any  new  pioneering  field,  the  dissemination  of  the  work  is  of  key  importanee.  Besides 
the  aforementioned  eonferenees,  the  Pis  have  aetively  presented  their  works  in  different  venues. 
PI  Jornet  has  taught  three  invited  eourses  and  three  invited  seminars  on  THz  eommunieations  in 
national  and  international  universities.  In  addition,  he  organized  and  hosted  the  first  Joint  Aeademia 
and  Industry  Panel  on  THz  Communieations,  at  IEEE  GLOBECOM  2015,  one  of  the  premiere 
eonferenees  of  the  IEEE  Communieations  Soeiety.  Co-PI  Einarsson  presented  results  at  the  2016 
IRMMW-THz  eonferenee  and  at  an  invited  seminar  at  Binghamton  University,  and  his  student 
supported  by  this  projeet  presented  her  work  at  the  2016  Ameriean  Physieal  Soeiety  Mareh  Meeting. 
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10  LIST  OF  ACRONYMS 


AFOSR 

Air  Eorce  Office  of  Scientific  Research 

AFRL 

Air  Eorce  Research  Eaboratory 

BER 

Bit  Error  Rate 

CSMA 

Carrier-Sense  Multiple  Access 

CTMA 

Continuous-time  Moving  Average 

FER 

Erame  Error  Rate 

MAC 

Medium  Access  Control 

PDE 

Probability  Density  Eunction 

PSD 

Power  Spectral  Density 

SER 

Symbol  Error  Rate 

SNR 

Signal-to-Noise  Ratio 

SPP 

Surface  Plasmon  Polariton 

USAE 

United  States  Air  Eorce 

VCD 

Voltage-controlled  Delay 
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